forum.alglib.net

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')

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