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 systemCode:
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;
}