forum.alglib.net :: View topic - C++ linlsqrsolvesparse() not working as expected.

#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;
}

Author:	DLoscos [ Fri Mar 01, 2019 3:21 pm ]
Post subject:	C++ linlsqrsolvesparse() not working as expected.
I am trying to use the linlsqrsolvesparse() to solve a non symmetric 25x25 linear system, but I keep getting termination code 4: \|\|A^TRk\|\|/(\|\|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_SIZEMESH_SIZE #define MESH_SIZE2 (2MESH_SIZE)-1 #define MESH_WIDTH (PLATE_WIDTH/(MESH_SIZE-1)) #define Origq 20000000 #define q -(OrigqMESH_WIDTHMESH_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 = 2MESH_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= 2MESH_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]= (KmV[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]= (KmV[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= 2MESH_SIZE/5; double tot= (MESH_SIZE2-1)(MESH_SIZE-1); double disp = maxCtot; 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 = (iMESH_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 = (iMESH_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, (iMESH_SIZE)+MESH_SIZE-1,(iMESH_SIZE)+MESH_SIZE-1, 1.0); bb[(iMESH_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 = (iMESH_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 = (iMESH_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 = (iMESH_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 = (iMESH_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 = (iMESH_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 = (iMESH_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; }

forum.alglib.net http://forum.alglib.net/

C++ linlsqrsolvesparse() not working as expected. http://forum.alglib.net/viewtopic.php?f=2&t=3947	Page 1 of 1

Page 1 of 1	All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/