00001
00002
00004
00005 #include "matrice.h"
00006 #include "check.h"
00007
00008 #define ABS(X) (X>0 ? X : -X)
00009
00010
00012
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
00063
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
00076
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
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
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
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
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
00362 }
00363
00365
00366
00367
00368
00369
00370 LUdek::LUdek(void)
00371 {
00372
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
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
00394 static int i,j,dim,poc=1;
00395 static Matrica R(6,6), RR(6,6), I(6,6), c(6,1);
00396
00397 CHECK((dim=M.GetCol()) == M.GetRow());
00398
00399
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
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
00435
00436
00437
00438
00439
00440
00441 for(i=0;i<n-1;i++)
00442 { if(pivot==1)
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 {
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
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
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 {
00479 b.data[P[i]][0] /= A.data[P[i]][i];
00480 for(j=0;j<i;j++)
00481
00482 b.data[P[j]][0] -= A.data[P[j]][i]*b.data[P[i]][0];
00483 }
00484 return 1;
00485 }