• Main Page
  • Modules
  • Classes
  • Files
  • File List

D:/Projekt/ECF_trunk/examples/iprojekt/matrice.cpp

00001 
00002 // Implementacija razreda Matrica, SLAJ, LUDek
00004 
00005 #include "matrice.h"
00006 #include "check.h"
00007 
00008 #define ABS(X)  (X>0 ? X : -X)
00009 //#include <axh.h>  // za assert
00010 
00012 // razred Matrica
00013 
00014 void Matrica::Free(void)
00015 {
00016     static int i;
00017     if(data==NULL)
00018         return;
00019     for(i=0;i<row;i++)
00020         delete[] data[i];
00021     delete[] data;
00022 }
00023 
00024 void Matrica::Take(void)
00025 {
00026     static int i;
00027     CHECK(row>0 && col>0);
00028     data=new double*[row];
00029     for(i=0;i<row;i++)
00030         data[i]=new double[col];
00031 }
00032 
00033 Matrica::Matrica(void)
00034 {
00035     col=row=0;
00036     data=NULL;
00037 }
00038 
00039 Matrica::Matrica(int _row,int _col=1)
00040 {
00041     col=_col;
00042     row=_row;
00043     Take();
00044 }
00045 
00046 void Matrica::Init(int _row, int _col)
00047 {
00048     Free();
00049     col=_col;
00050     row=_row;
00051     Take();
00052 }
00053 
00054 Matrica::Matrica(Matrica &source)
00055 {
00056     static int i;
00057     row=source.row;
00058     col=source.col;
00059     Take();
00060     for(i=0;i<row;i++)
00061         memcpy(data[i],source.data[i],col*sizeof(double));
00062         /*for(int j=0;j<col;j++)
00063             data[i][j]=source.data[i][j];*/
00064 }
00065 
00066 Matrica& Matrica::operator =(Matrica& source)
00067 {
00068     static int i;
00069     if(data!=NULL) Free();
00070     row=source.row;
00071     col=source.col;
00072     Take();
00073     for(i=0;i<row;i++)
00074         memcpy(data[i],source.data[i],col*sizeof(double));
00075         /*for(int j=0;j<col;j++)
00076             data[i][j]=source.data[i][j];*/
00077     return *this;
00078 }
00079 
00080 Matrica::~Matrica()
00081 {
00082     Free();
00083 }
00084 
00085 Matrica Matrica::operator +(Matrica &op2)
00086 {
00087     Matrica R(row,col);
00088     CHECK(row==op2.row && col==op2.col && row>0);
00089     for(int i=0;i<row;i++)
00090         for(int j=0;j<col;j++)
00091             R.data[i][j]=data[i][j]+op2.data[i][j];
00092     return R;
00093 }
00094 
00095 void Matrica::operator +=(Matrica &op2)
00096 {
00097     CHECK(row==op2.row && col==op2.col && row>0);
00098     for(int i=0;i<row;i++)
00099         for(int j=0;j<col;j++)
00100             data[i][j] += op2.data[i][j];
00101 }
00102 
00103 void Matrica::operator -=(Matrica &op2)
00104 {
00105     CHECK(row==op2.row && col==op2.col && row>0);
00106     for(int i=0;i<row;i++)
00107         for(int j=0;j<col;j++)
00108             data[i][j] -= op2.data[i][j];
00109 }
00110 
00111 Matrica Matrica::operator -(const Matrica &op2)
00112 {
00113     Matrica R(row,col);
00114     CHECK(row==op2.row && col==op2.col && row>0);
00115     for(int i=0;i<row;i++)
00116         for(int j=0;j<col;j++)
00117             R.data[i][j]=data[i][j]-op2.data[i][j];
00118     return R;
00119 }
00120 
00121 Matrica Matrica::operator *(const Matrica &op2)
00122 {
00123     Matrica R(row,op2.col);
00124     CHECK(col==op2.row);
00125     for(int i=0;i<row;i++)
00126         for(int j=0;j<op2.col;j++)
00127         {   R.data[i][j]=0;
00128             for(int k=0;k<col;k++)
00129                 R.data[i][j]+=data[i][k]*op2.data[k][j];
00130         }
00131     return R;
00132 }
00133 
00134 Matrica Matrica::operator *(const double &scalar)
00135 {
00136     Matrica R(row,col);
00137     for(int i=0;i<row;i++)
00138         for(int j=0;j<col;j++)
00139             R.data[i][j]=this->data[i][j]*scalar;
00140     return R;
00141 }
00142 
00143 // translacija
00144 Matrica Matrica::operator ~(void)
00145 {
00146     CHECK(data!=NULL);
00147     Matrica P(col, row);
00148     int i,j;
00149     for(i=0;i<row;i++)
00150         for(j=0;j<col;j++)
00151             P.data[j][i]=data[i][j];
00152     return P;
00153 }
00154 
00155 Matrica Matrica::operator /(const double &scalar)
00156 {
00157     Matrica R(row,col);
00158     int i,j;
00159 
00160     for (i=0;i<row; i++)
00161         for (j=0;j<col; j++)
00162             R.Set(i,j,data[i][j]/scalar);
00163     return R;
00164 }
00165 
00166 // transponira matricu: ako je parametar 0, transponira samo kopiju, inace sebe samu
00167 Matrica Matrica::Transpose(int transpose)
00168 {
00169     CHECK(data!=NULL);
00170     Matrica P(col, row);
00171     int i,j;
00172     for(i=0;i<row;i++)
00173         for(j=0;j<col;j++)
00174             P.data[j][i]=data[i][j];
00175     if(!transpose)
00176         return P;
00177     // ako i this treba transponirati
00178     i=row; row=col; col=i;
00179     Free(); Take();
00180     for(i=0;i<row;i++)
00181         for(j=0;j<col;j++)
00182             data[i][j]=P.data[i][j];
00183     return (*this);
00184 }
00185 
00186 double Matrica::Get(int _row,int _col) const
00187 {
00188     CHECK(data!=NULL && _row>=0 && _row<row && _col>=0 && _col<col);
00189     return data[_row][_col];
00190 }
00191 
00192 double Matrica::Get(int _row) const
00193 {
00194     CHECK(data!=NULL && _row>=0 && _row<row);
00195     return data[_row][0];
00196 }
00197 
00198 double Matrica::Set(int _row,int _col, double value)
00199 {
00200     CHECK(data!=NULL && _row>=0 && _row<=row && _col>=0 && _col<=col);
00201     data[_row][_col]=value;
00202     return value;
00203 }
00204 
00205 double Matrica::Set(int _row,double value)
00206 {
00207     CHECK(data!=NULL && _row>=0 && _row<=row);
00208     data[_row][0]=value;
00209     return value;
00210 }
00211 
00212 int Matrica::Load(const char* fname)
00213 {
00214     double value;
00215     FILE *fp;
00216     CHECK((fp=fopen(fname,"r"))!=NULL);
00217     if(!(fscanf(fp,"%d %d",&row,&col)!=EOF))
00218         exit(1);
00219     if(data!=NULL)
00220         Free();
00221     Take();
00222     for(int i=0;i<row*col;i++)
00223     {   if(!(fscanf(fp," %lf ",&value)!=EOF))
00224             exit(1);
00225         Set(i/col,i%col,value);
00226     }
00227     fclose(fp);
00228     return 1;
00229 }
00230 
00231 void Matrica::Save(const char* fname)
00232 {
00233     FILE *fp;
00234     if((fp=fopen(fname,"w+"))==NULL)
00235         exit(1);
00236     fprintf(fp,"%d %d\n",row,col);
00237     for(int i=0;i<row;i++)
00238     {   for(int j=0;j<col;j++)
00239             fprintf(fp," %g ",data[i][j]);
00240         fprintf(fp,"\n");
00241     }
00242     fclose(fp);
00243 }
00244 
00245 void Matrica::SaveLong(const char* fname)
00246 {
00247     FILE *fp;
00248     if((fp=fopen(fname,"w+"))==NULL)
00249         exit(1);
00250     fprintf(fp,"%d %d\n",row,col);
00251     for(int i=0;i<row;i++)
00252     {   for(int j=0;j<col;j++)
00253             fprintf(fp," %lf ",data[i][j]);
00254         fprintf(fp,"\n");
00255     }
00256     fclose(fp);
00257 }
00258 
00259 void Matrica::Show(void)
00260 {
00261     int i,j;
00262     printf("Rows: %d Cols: %d\n",row,col);
00263     for(i=0;i<row;i++)
00264     {   for(j=0;j<col;j++)
00265             printf("%G\t",data[i][j]);
00266         printf("\n");
00267     }
00268 }
00269 
00270 void Matrica::Reset(int _row,int _col)
00271 {
00272     CHECK(_row>0 && _col>0);
00273     Free();
00274     row=_row;
00275     col=_col;
00276     Take();
00277     for(int i=0; i<row; i++)
00278         for(int j=0; j<col; j++)
00279             this->data[i][j] = 0;
00280 }
00281 
00282 int Matrica::CopyColumn(int colnum,Matrica &V)
00283 {
00284     int i;
00285     CHECK(colnum>-1 && colnum<col);
00286     V.Reset(row,1);
00287     for(i=0;i<row;i++)
00288         V.Set(i,data[i][colnum]);
00289     return 1;
00290 }
00291 
00292 int Matrica::SetColumn(int colnum,Matrica &v)
00293 {
00294     int i;
00295     CHECK(v.GetRow()==row);
00296     CHECK(colnum>-1 && colnum<col);
00297     for(i=0;i<row;i++)
00298         data[i][colnum]=v.Get(i);
00299     return 1;
00300 }
00301 
00302 double Matrica::Norm(void)
00303 {
00304     double r=0;
00305     for(int i=0;i<row;i++)
00306         r+= data[i][0]*data[i][0];
00307     return r;
00308 }
00309 
00310 Matrica Matrica::operator+ (double broj)
00311 {   Matrica R(*this);
00312     for(int i=0; i<row; i++)
00313         for(int j=0; j<col; j++)
00314             R[i][j] += broj;
00315     return R;
00316 }
00317 
00318 
00320 // razred SLAJ
00321 
00322 void SLAJ::Load_A(char* fname)
00323 {
00324     A.Load(fname);
00325 }
00326 
00327 void SLAJ::Load_b(char* fname)
00328 {
00329     b.Load(fname);
00330 }
00331 
00332 void SLAJ::Save_A(char* fname)
00333 {
00334     A.Save(fname);
00335 }
00336 
00337 void SLAJ::Save_b(char* fname)
00338 {
00339     b.Save(fname);
00340 }
00341 
00342 void SLAJ::Set_A(Matrica &_A)
00343 {
00344     A.Reset(_A.GetRow(),_A.GetCol());
00345     A=_A;
00346 }
00347 
00348 void SLAJ::Set_b(Matrica &_b)
00349 {
00350     CHECK(_b.GetCol()==1);
00351     b=_b;
00352 }
00353 
00354 void SLAJ::GetSolution(Matrica &X)
00355 {
00356     static int row,i;
00357     row=b.GetRow();
00358     X.Reset(row,1);
00359     for(i=0;i<row;i++)
00360         X.data[i][0]=b.data[i][0];
00361     //X=b;
00362 }
00363 
00365 // razred LUdek
00366 
00367 // za potrebe DASGUPTA6 promijenio sam sve oznaceno sa !!!
00368 
00369 
00370 LUdek::LUdek(void)
00371 {
00372 //  P=NULL; !!!
00373     P=new int[9];
00374     for(int i=0;i<9;i++)
00375         P[i]=i;
00376 }
00377 
00378 int LUdek::Solve(Matrica &A,Matrica &X,Matrica &b)
00379 {
00380     Set_A(A);
00381     Set_b(b);
00382     Decompose();
00383     ForwardSubst();
00384     BackwardSubst();
00385     //GetSolution(X);
00386     for(int i=0;i<n;i++)
00387         X.data[i][0]=b.data[P[i]][0];
00388     return 1;
00389 }
00390 
00391 Matrica LUdek::Invert(Matrica &M)
00392 {
00393     // rjesavamo sustav M*M!=I
00394     static int i,j,dim,poc=1;
00395     static Matrica R(6,6), RR(6,6), I(6,6), c(6,1);
00396     //static Matrica R(3,3), RR(3,3), I(3,3), c(3,1);
00397     CHECK((dim=M.GetCol()) == M.GetRow());
00398 //  Matrica R(dim,dim), RR(dim,dim), I(dim,dim), c(dim,1); !!!
00399     // I jedinicna
00400     if(poc==1)
00401     {   for(i=0;i<dim;i++)
00402         {   for(j=0;j<dim;j++)
00403                 I.Set(i,j,0);
00404             I.Set(i,i,1);
00405         }
00406         poc=0;
00407     }
00408 
00409     Set_A(M);
00410     Decompose();
00411 
00412     for(i=0;i<dim;i++)
00413     {   I.CopyColumn(i,c);
00414         Set_b(c);
00415         ForwardSubst();
00416         BackwardSubst();
00417         GetSolution(c);
00418         R.SetColumn(i,c);
00419     }
00420 
00421     //vratimo retke ako je bilo pivotiranja
00422     for(i=0;i<n;i++)
00423         for(j=0;j<n;j++)
00424             RR.data[i][j]=R.data[P[i]][j];
00425 
00426     return RR;
00427 }
00428 
00429 
00430 int LUdek::Decompose(void)
00431 {
00432     static int i,j,k,p;
00433     n=A.GetRow();
00434 /*  if(P!=NULL) !!!
00435         delete[] P;
00436     P=new int[n];
00437     for(i=0;i<n;i++)
00438         P[i]=i;
00439 */
00440     
00441     for(i=0;i<n-1;i++)
00442     {   if(pivot==1)        // pivotiranje
00443         {   p=i;
00444             for(j=i+1;j<n;j++)
00445                 if(ABS(A.Get(P[j],i))>ABS(A.Get(P[p],i)))
00446                     p=j;
00447             if(A.Get(P[p],i)==0)
00448                 CHECK(0);
00449             k=P[i];
00450             P[i]=P[p];
00451             P[p]=k;
00452         }
00453         for(j=i+1;j<n;j++)
00454         {   //A.Set(P[j],i,A.Get(P[j],i)/A.Get(P[i],i));
00455             A.data[P[j]][i]=A.data[P[j]][i]/A.data[P[i]][i];
00456             for(k=i+1;k<n;k++)
00457                 //A.Set(P[j],k,A.Get(P[j],k)-A.Get(P[j],i)*A.Get(P[i],k));
00458                 A.data[P[j]][k]-= A.data[P[j]][i]*A.data[P[i]][k];
00459         }
00460     }
00461     return 1;
00462 }
00463 
00464 int LUdek::ForwardSubst(void)
00465 {
00466     static int i,j;
00467     for(i=0;i<n-1;i++)
00468         for(j=i+1;j<n;j++)
00469             //b.Set(P[j],0,b.Get(P[j],0)-A.Get(P[j],i)*b.Get(P[i],0));
00470             b.data[P[j]][0]-= A.data[P[j]][i]*b.data[P[i]][0];
00471     return 1;
00472 }
00473 
00474 int LUdek::BackwardSubst(void)
00475 {
00476     static int i,j;
00477     for(i=n-1;i>=0;i--)
00478     {   //b.Set(P[i],0,b.Get(P[i],0)/A.Get(P[i],i));
00479         b.data[P[i]][0] /= A.data[P[i]][i];
00480         for(j=0;j<i;j++)
00481             //b.Set(P[j],0,b.Get(P[j],0)-A.Get(P[j],i)*b.Get(P[i],0));
00482             b.data[P[j]][0] -= A.data[P[j]][i]*b.data[P[i]][0];
00483     }
00484     return 1;
00485 }

Generated on Thu Jul 10 2014 14:13:41 for ECF by  doxygen 1.7.1