math - LU decomposition in c++ not working on big matrices -

update: here solution, added scalar each row bring underflow, overflow under control. guys.

i have been working on lu decomposition in c++ 1 day decompose , solve large sparse matrix. found code , modified own use not work on big matrices. works on matrices size 5 5. need work matrices of size 100 100 or more. have checked solutions in mat-lab , code gives totally wrong results. feel problem comes division in code, , if so, suggestions how solve appreciated ans appreciated.

here code.

updated:

 #include <algorithm>  // **  * end ***  /*  * ludecomp.cpp  #include <stdlib.h>  #include <stdio.h>  #include <math.h>  #include <iostream>  #include <fstream>  #include <string.h>  #include <iomanip>   #include "ludecomp.h"   using namespace std;   ludecomp::ludecomp()  {  }   void ludecomp::h_pivot_decomp(int mat1, double a[], int p[], int q[])  { int = 0, j = 0, k = 0; int n = mat1; int pi = 0, pj = 0, tmp = 0; double max = 0.0; double ftmp = 0.0; //stores scaling of each row or column. double* vv = new double[mat1 + 5];  //loop on rows toget implicit scaling information. max = 0.0; (i = 0; < n; i++) {     (j = 0; j < n; j++)     {          if ((ftmp = fabs(a(i,j))) > max)         {             max=ftmp;         }     }     //no nonzero largest element.                 if (max == 0.0)                 {                     throw("singular matrix in ludcmp");                 }                 //save scaling.                 vv[i]=1.0/max;              }          // k element determines pivot element in thereby         // determining submatrix starting @ upper left corner of matrix. (k = 0; k < n; k++) {      // pi: stores row needing swapped.     // pj: stores column needing swapped.     // max: makes 0 element in matrix tiny number.     pi = -1, pj = -1, max = tiny;      //find pivot in submatrix a(k:n,k:n) finding absolute value of biggest element.     (i = k; < n; i++)     {         (j = k; j < n; j++)         {             //j = k;             ftmp = vv[i] * fabs(a(i,j));             // decides if current max bigger current element.                     if (ftmp>max)                     {                         max = ftmp;                         // index of row being swapped.                         pi=i;                         // index of column being swapped.                         pj=j;                     }                 }             }      {         // stores permutation of row swaps.         tmp = p[k];         p[k] = p[pi];         p[pi] = tmp;     }      //swaps scalling factor if needed.     if (k != pi)     {         vv[pi] = vv[k];         cout << "scaling factor: " << vv[pi] << endl;     }      // swaps indicated rows move max pivot     // element of submatrix k place.     (j = 0; j < n; j++)     {         // k , pi index stays same row         // number stays same, j changes iterate threw row.         ftmp = a(k,j);         a(k,j)=a(pi,j);         a(pi,j)=ftmp;         //cout << a(k,j) << " , " << a(pi,j) << endl;             }      {         // stores permutation of column swaps.         tmp = q[k];         q[k] = q[pj];         q[pj] = tmp;         //cout << q[k] << " , " << q[pj] << endl;     }      // swaps indicated columns move max pivot     // element of submatrix k place.     (i = 0; < n; i++)     {         // k , pj index stays same column         // number stays same, changes iterate threw column.         ftmp = a(i,k);         a(i,k)=a(i,pj);         a(i,pj)=ftmp;         //cout << a(i,k) << " , " << a(i,pj) << endl;             }         // end pivot     cout << fixed << showpoint;     cout << setprecision(20);     // check pivot size , decompose     if ((fabs(a(k,k))>tiny))     {         (i=k+1;i<n;i++)         {             // column normalisation, first element under pivot k row i.              ftmp=a(i,k)/=a(k,k);              cout << "k,k " <<a(k,k) << " , " << endl;             // rest of row i.             (j=k+1;j<n;j++)             {                 //a(ik)*a(kj) subtracted lower right submatrix elements                 a(i,j)-=(ftmp*a(k,j));                 //cout <<"i,j "<< a(i,j) << endl;             }         }     }  }     //end decompose (i = 0; < n; i++) {     (j = 0; j < n; j++)     {          cout << a(i,j)<<" ";     }     cout << endl; }  }   void ludecomp::h_solve(int mat1, double a[], double x[], int p[], int q[])  { // forward substitution; see  golub, van loan 96 // , see http://www.cs.rutgers.edu/~richter/cs510/completepivoting.pdf int = 0, ii = 0, j = 0; double ftmp = 0.0; double* xtmp = new double[mat1 + 5];  cout << fixed << showpoint; cout << setprecision(4);  // swap rows // put vector should using permutations row swapping. (i = 0; < mat1; i++) {     xtmp[i] = x[p[i]]; //value should here     //cout << xtmp[i] << endl; }  // ly=b (i = 0; < mat1; i++) {     ftmp = xtmp[i];     if (ii != 0)         (j = ii - 1; j < i; j++)             ftmp -= a(i,j)*xtmp[j];              else if (ftmp!=0.0)             ii=i+1;      xtmp[i] = ftmp;     //cout << xtmp[i] << endl; }  // backward substitution // partially taken sourcebook on parallel computing p577 // solves ux=y cout << "xtmp " << xtmp[mat1 - 1] << " " << a(mat1-1,mat1-1)<< endl; xtmp[mat1 - 1] /= a(mat1-1,mat1-1); //cout << xtmp[mat1 - 1] << endl; (i = mat1 - 2; >= 0; i--) {     ftmp = xtmp[i];     //cout << "ftmp " << ftmp << endl;     (j = + 1; j < mat1; j++)     {         ftmp -= a(i,j)*xtmp[j];         //cout << "ftmp in "<<ftmp << endl;     }      xtmp[i] = (ftmp) / a(i,i);  }      // last bit     // swap columns     // takes final answer , puts proper order     // using permutations column swapping. (i = 0; < mat1; i++) {     x[q[i]] = xtmp[i]; }  delete xtmp;  }   // method output lu decomposition.  void ludecomp::output(unsigned int mat1, double a[], double b[])  {  // pivot array's permutation vectors. int* p_pivot = new int[mat1 + 5]; int* q_pivot = new int[mat1 + 5];  // sets elements in permutation vectors receive permutations. // p_pivot row permutations , initialized {0,1,...,r}; // q_pivot column permutations , initialized {0,1,...,r}; (unsigned int = 0; < mat1; i++) {     p_pivot[i] = i;     q_pivot[i] = i; }  // call decomposition method passing (size,matrix decomposed, not used,   not used). h_pivot_decomp(mat1, a, p_pivot, q_pivot);  // call solve passing (size, matrix in lu form, b vector, not used, not used). h_solve(mat1, a, b, p_pivot, q_pivot);  // have solution. // used file output. ofstream outfile; // allow appenending file created. outfile.open("outsolmatrix.txt");  // sets precision of output file. outfile << fixed << showpoint; outfile << setprecision(4);  // output results file answer {0,1,...,n}. (unsigned int = 0; < mat1; i++) {     outfile << << " " << b[i] << endl; }  outfile << "end" << endl;  delete p_pivot; delete q_pivot;  outfile.close();  } 

the h file here if need see it.

#ifndef ludecomp_h_ #define ludecomp_h_  class ludecomp {  public:  #define a(i,j) a[(i)*mat1+(j)]  const static double tiny = 1e-20;  ludecomp();  void h_pivot_decomp(int mat1, float *a, int *p, int *q); void h_solve(int mat1, float *a, float *x, int *p, int *q); void output(unsigned int mat1, float *a, float *b);  private:  };  #endif /* ludecomp_h_ */ 

thanks again, let me know if guys need see else.

try replace float double. float not type complex computations, because small (only 4 bytes), not exact (only 7-8 digits). better use double (8 bytes, 15-16 digits). if need more exact computatuions, have use data structures values (like in matlab, maple , other systems). see arbitrary-precision arithmetic.