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 |
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^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; } |
Page 1 of 1 | All times are UTC |
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/ |