forum.alglib.net

ALGLIB forum
It is currently Sun Oct 20, 2019 2:14 pm

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 1 post ] 
Author Message
 Post subject: C++ linlsqrsolvesparse() not working as expected.
PostPosted: Fri Mar 01, 2019 3:21 pm 
Offline

Joined: Fri Mar 01, 2019 3:04 pm
Posts: 1
I am trying to use the linlsqrsolvesparse() to solve a non symmetric 25x25 linear system, but I keep getting termination code 4: ||A^T*Rk||/(||A||*||Rk||)<=EpsA
With a null vector as a solution.

I have tested the same system in matlab and it was able to compute it with no problems, so I don't think my problem is ill conditioned.

To check if my code was written correctly I solved the example from the linlgcolvesparse() with the GC method and then with the LSQR method. The results matched and everything worked properly.

I have throughtly read the API and played with different conditioning and parameters for the solver but I always get the 0 solution with termination code 4.
Any ideas on why is this happening or how to fix it?

matlab solution for the system
Code:
A =[-2 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 -3 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 -3 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 -3 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 -3 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 -3 1 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 1 -3 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 -3 1 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -4 1 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -3 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 -2 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -3 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -3 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -3 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 -2]

B=[-31.2500,-31.2500,-31.2500,-31.2500,293.0000,-31.2500,-31.2500,-31.2500,-31.2500,293.0000,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500,-31.2500]
T = reshape(A\B',5,5)
imagesc(T')


main.cpp for reference:
Code:
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "solvers.h"
#include <iostream>
#include <string>

#define MESH_SIZE             5
#define maxC                  0.4
#define PLATE_WIDTH                0.005
#define SINK_WIDTH             0.002
#define Ts                   293
#define Km                   65.0
#define Kp                  0.2
#define TEMP_SIZE               MESH_SIZE*MESH_SIZE
#define MESH_SIZE2            (2*MESH_SIZE)-1
#define MESH_WIDTH             (PLATE_WIDTH/(MESH_SIZE-1))
#define Origq             20000000
#define q                 -(Origq*MESH_WIDTH*MESH_WIDTH)
#define K()               1


using namespace std;
using namespace alglib;

void initV(double V[MESH_SIZE2][MESH_SIZE],double fill);
void computeK(double V[MESH_SIZE2][MESH_SIZE],double K[MESH_SIZE2][MESH_SIZE]);
double Vcon(double V[MESH_SIZE2][MESH_SIZE]);
double V0();
real_1d_array generateSystem(double K[MESH_SIZE2][MESH_SIZE], int metal, sparsematrix M);
string array2string(double bb[TEMP_SIZE]);

int main(int argc, char *argv[]) {
   
   double fill = V0();
  double VMesh [MESH_SIZE2][MESH_SIZE] = {0};
  double KMesh [MESH_SIZE2][MESH_SIZE] = {0};
   int metal = 2*MESH_SIZE/5;
  initV(VMesh,fill);
   computeK(VMesh,KMesh);


   //TEST SYMMETRIC SYSTEM FROM API
   sparsematrix a;
    sparsecreate(5, 5, a);
    sparseset(a, 0, 0, 5.0);
    sparseset(a, 0, 1, 1.0);
    sparseset(a, 1, 0, 1.0);
    sparseset(a, 1, 1, 7.0);
    sparseset(a, 1, 2, 2.0);
    sparseset(a, 2, 1, 2.0);
    sparseset(a, 2, 2, 8.0);
    sparseset(a, 2, 3, 1.0);
    sparseset(a, 3, 2, 1.0);
    sparseset(a, 3, 3, 4.0);
    sparseset(a, 3, 4, 1.0);
    sparseset(a, 4, 3, 1.0);
    sparseset(a, 4, 4, 4.0);

    sparseconverttocrs(a);

    real_1d_array b = "[7,17,14,10,6]";
    lincgstate s;
    lincgreport rep;
    real_1d_array x;
    lincgcreate(5, s);
    lincgsolvesparse(s, a, true, b);
    lincgresults(s, x, rep);

    printf("%d\n", int(rep.terminationtype)); // EXPECTED: 1
      printf("%d\n", int(rep.iterationscount));
    printf("%s\n", x.tostring(2).c_str()); // EXPECTED: [1.000,2.000,1.000,2.000,1.000]

   //SAME TEST SYSTEM WITH   LSQR SOLVER
   linlsqrstate solver2;
  linlsqrreport report2;
  real_1d_array T2;
   linlsqrcreate(5,5, solver2);
  linlsqrsolvesparse(solver2, a, b);
   linlsqrresults(solver2, T2, report2);

  printf("%d\n", int(report2.terminationtype));
   printf("%d\n", int(report2.iterationscount));
   printf("%d\n", int(report2.nmv));
  printf("%s\n", T2.tostring(4).c_str());


   //ACTUAL NON-SYMMETRIC SYSTEM WITH LSQR SOLVER
   sparsematrix M;
   sparsecreate(TEMP_SIZE, TEMP_SIZE, M);
  real_1d_array Q = generateSystem(KMesh,metal,M);
   printf("%s\n", Q.tostring(4).c_str());
   sparseconverttocrs(M);
   linlsqrstate solver;
//linlsqrsetcond(solver,0.000000000000000001,0,0);
//linlsqrsetprecunit(solver);
  linlsqrreport report;
  real_1d_array T;
   linlsqrcreate(TEMP_SIZE,TEMP_SIZE, solver);
  linlsqrsolvesparse(solver, M, Q);
   linlsqrresults(solver, T, report);

  printf("%d\n", int(report.terminationtype)); // EXPECTED: 4
   printf("%d\n", int(report.iterationscount)); // GETTING: 0
   printf("%d\n", int(report.nmv));
  printf("%s\n", T.tostring(4).c_str());

  return 0 ;
}


void initV(double V[MESH_SIZE2][MESH_SIZE],double fill){
   int metal= 2*MESH_SIZE/5;
   cout << "metal = " << metal << endl;
   int i = 0;
   if (metal%2 == 0){
      while (i<metal){
         for (int j=0; j<MESH_SIZE-2; j++){
            V[i][j] = fill;
         }
         V[i][MESH_SIZE-2] = 1;
         i++;
         for (int j=0; j<MESH_SIZE-1; j++){
            V[i][j] = fill;
         }
         V[i][MESH_SIZE-1] = 1;
         i++;
      }
      while (i<MESH_SIZE2-2){
         for (int j=0; j<MESH_SIZE-1; j++){
            V[i][j] = fill;
         }
         i++;
         for (int j=0; j<MESH_SIZE; j++){
            V[i][j] = fill;
         }
         i++;
      }
      for (int j=0; j<MESH_SIZE-1; j++){
            V[i][j] = fill;
      }
   }else{
         for (int j=0; j<MESH_SIZE-2; j++){
                V[i][j] = fill;
            }
         V[i][MESH_SIZE-2] = 1;
         i++;
      while (i<metal-1){
         for (int j=0; j<MESH_SIZE-1; j++){
            V[i][j] = fill;
         }
         V[i][MESH_SIZE-1] = 1;
         i++;
         for (int j=0; j<MESH_SIZE-2; j++){
            V[i][j] = fill;
         }
         V[i][MESH_SIZE-2] = 1;
         i++;
      }
      while (i<MESH_SIZE2-1){
         for (int j=0; j<MESH_SIZE; j++){
            V[i][j] = fill;
         }
         i++;
         for (int j=0; j<MESH_SIZE-1; j++){
            V[i][j] = fill;
         }
         i++;
      }
   }
}

void computeK(double V[MESH_SIZE2][MESH_SIZE],double K[MESH_SIZE2][MESH_SIZE]){
for (int i=0; i<MESH_SIZE2; i=i+2){
      for (int j=0; j<MESH_SIZE-1; j++){
        K[i][j]= (Km*V[i][j]) + Kp*(1-V[i][j]);
      }
   }
   for (int i=1; i<MESH_SIZE2; i=i+2){
      for (int j=0; j<MESH_SIZE; j++){
        K[i][j]= (Km*V[i][j]) + Kp*(1-V[i][j]);
      }
   }
}

double Vcon(double V[MESH_SIZE2][MESH_SIZE]){
   double con = 0;
  double tot = (MESH_SIZE2-1)*(MESH_SIZE-1);
   for (int j=0; j<MESH_SIZE-1; j++){
        con += (V[0][j]+V[MESH_SIZE2-1][j])/2;
   }
   for (int i=2; i<MESH_SIZE2-2; i=i+2){
      for (int j=0; j<MESH_SIZE-1; j++){
        con += V[i][j];
      }
   }
   for (int i=1; i<MESH_SIZE2-1; i=i+2){
      con += (V[i][0] + V[i][MESH_SIZE-1])/2;
      for (int j=1; j<MESH_SIZE-1; j++){
        con += V[i][j];
      }
   }
  return con/tot;
}

double V0(){
   int metal= 2*MESH_SIZE/5;
   double tot= (MESH_SIZE2-1)*(MESH_SIZE-1);
  double disp = maxC*tot;
   disp = disp-(metal/2.0);
   return disp/(tot-(metal/2));
}

real_1d_array generateSystem(double K[MESH_SIZE2][MESH_SIZE], int metal, sparsematrix M){
    sparsecreate(TEMP_SIZE, TEMP_SIZE, M);
      double bb[TEMP_SIZE] = {0};
      double qvalue =  q;
      int t_index;
      //Top border
      sparseset(M,0,0, -(K()+K()));//corner
      sparseset(M,0,1, K());//right
      sparseset(M,0,MESH_SIZE, K());//down
      bb[0] = q;
      for (int j = 1; j < MESH_SIZE-1; j++){
                  sparseset(M,j,j, -(K()+K()+K()));
                  sparseset(M,j,j-1, K());//left
                  sparseset(M,j,j+1, K());//right
                  sparseset(M,j,j+MESH_SIZE, K());//down
                  bb[j] = qvalue;
      }
      sparseset(M,MESH_SIZE-1,MESH_SIZE-1, 1.0);//sink corner
      bb[MESH_SIZE-1] = Ts;
   
      //equation for the rows with heatsink points   
      for(int i = 1; i < metal; i++){
            //Left Border
                  t_index = (i*MESH_SIZE);
                  sparseset(M,t_index,t_index, -(K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index+1, K());//right
                  sparseset(M,t_index,t_index+MESH_SIZE, K());//down
                  bb[t_index] = qvalue;
            //Interior points
            for (int j = 1; j < MESH_SIZE-1; j++){
                  t_index = (i*MESH_SIZE)+j;
                  sparseset(M,t_index,t_index, -(K()+K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index-1, K());//left
                  sparseset(M,t_index,t_index+1, K());//right
                  sparseset(M,t_index,t_index+MESH_SIZE, K());//down
                  bb[t_index] = qvalue;
            }
            //heat sink points
            sparseset(M, (i*MESH_SIZE)+MESH_SIZE-1,(i*MESH_SIZE)+MESH_SIZE-1, 1.0);
            bb[(i*MESH_SIZE)+MESH_SIZE-1] = Ts;
      }
      //equation for the rows without heatsink points   
      for(int i = metal; i < MESH_SIZE -1; i++){
            //Left Border
                  t_index = (i*MESH_SIZE);
                  sparseset(M,t_index,t_index, -(K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index+1, K());//right
                  sparseset(M,t_index,t_index+MESH_SIZE, K());//down
                  bb[t_index] = qvalue;
            //Interior points
            for (int j = 1; j < MESH_SIZE-1; j++){
                  t_index = (i*MESH_SIZE)+j;
                  sparseset(M,t_index,t_index, -(K()+K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index-1, K());//left
                  sparseset(M,t_index,t_index+1, K());//right
                  sparseset(M,t_index,t_index+MESH_SIZE, K());//down
                  bb[t_index] = qvalue;
            }
            //right border points
                  t_index = (i*MESH_SIZE)+MESH_SIZE-1;
                  sparseset(M,t_index,t_index, -(K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index-1, K());//left
                  sparseset(M,t_index,t_index+MESH_SIZE, K());//down
                  bb[t_index] = qvalue;
      }
      //Bottom border
      int i = MESH_SIZE-1;
      //j = 0
      t_index = (i*MESH_SIZE);
      sparseset(M,t_index,t_index, -(K()+K()));//corner
      sparseset(M,t_index,t_index-MESH_SIZE, K());//up
      sparseset(M,t_index,t_index+1, K());//right
      bb[t_index] = qvalue;
      for (int j = 1; j < MESH_SIZE-1; j++){
                  t_index = (i*MESH_SIZE)+j;
                  sparseset(M,t_index,t_index, -(K()+K()+K()));
                  sparseset(M,t_index,t_index-MESH_SIZE, K());//up
                  sparseset(M,t_index,t_index-1, K());//left
                  sparseset(M,t_index,t_index+1, K());//right
                  bb[t_index] = qvalue;
      }
      int j = MESH_SIZE-1;
      t_index = (i*MESH_SIZE)+j;
      sparseset(M,t_index,t_index, -(K()+K()));//corner
      sparseset(M,t_index,t_index-MESH_SIZE, K());//up
      sparseset(M,t_index,t_index-1, K());//left
      cout << "Initial concentration matrix: " << endl;
   for (int i=0; i<TEMP_SIZE; i++){
      for (int j=0; j<TEMP_SIZE; j++){
        cout << sparseget(M,i,j) << " ";
      }
      cout <<  endl;
   }
      bb[t_index] = qvalue;

      //Actually build the q vector
      string s = array2string(bb);
      real_1d_array Qv = real_1d_array(s.c_str());
return Qv;
}


string array2string(double bb[TEMP_SIZE]){
   string s = "[";
   for (int i=0; i<TEMP_SIZE-1; i++){
      s = s + to_string(bb[i]) + ",";
   }
   s = s + to_string(bb[TEMP_SIZE-1]) + "]";
return s;
}


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 1 post ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 1 guest


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group