00001
00002
00003 #include "SparseMatrixF.h"
00004 #include "Array.h"
00005 #include "BigMatrix.h"
00006
00007 DSS_NAMESPASE_BEGIN
00008
00009 SparseMatrixF::SparseMatrixF()
00010 {
00011 a = NULL;
00012 ci = NULL;
00013 adr = NULL;
00014 neq = 0;
00015 Aoffset = 0;
00016 Coffset = 0;
00017 bJardaConvention = true;
00018 bLocalCopy = false;
00019 m_bIsSymmertric = false;
00020 m_eSparseOrientation = eCompressedColumns;
00021 }
00022
00023 SparseMatrixF::SparseMatrixF(unsigned long neq,double* a,unsigned long* ci,unsigned long *adr,ULONG Aofs ,ULONG Cofs ,bool JardaConvention ,bool bIsSymetric,eOrientation sparseOri)
00024 {
00025 this->a = a;
00026 this->neq = neq;
00027 this->adr = adr;
00028 this->ci = ci;
00029 Aoffset = Aofs;
00030 Coffset = Cofs;
00031 bJardaConvention = JardaConvention;
00032 bLocalCopy = false;
00033 m_bIsSymmertric = bIsSymetric;
00034 m_eSparseOrientation = sparseOri;
00035 }
00036
00037 SparseMatrixF::SparseMatrixF(IConectMatrix* )
00038 {
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 }
00060
00061 SparseMatrixF::~SparseMatrixF()
00062 {
00063 if (bLocalCopy)
00064 Delete();
00065 }
00066
00067 void SparseMatrixF::Delete()
00068 {
00069 if (a) {delete [] a;a=NULL;}
00070 if (ci) {delete [] ci;ci=NULL;}
00071 if (adr) {delete [] adr;adr = NULL;}
00072 neq = 0;
00073 bLocalCopy = false;
00074 }
00075
00076 void SparseMatrixF::Detach()
00077 {
00078 a=NULL;
00079 ci=NULL;
00080 adr = NULL;
00081 neq = 0;
00082 Aoffset = 0;
00083 Coffset = 0;
00084 }
00085
00086 void SparseMatrixF::CreateLocalCopy()
00087 {
00088 if (neq)
00089 {
00090 long neql = neq;
00091 if (bJardaConvention) neql++;
00092 long nnz = Nonzeros();
00093
00094 ULONG* old_adr = adr; adr = new ULONG[neql];
00095 Array::Copy((long*)old_adr,(long*)adr,neql);
00096
00097 ULONG* old_ci = ci; ci = new ULONG[nnz];
00098 Array::Copy((long*)old_ci,(long*)ci,nnz);
00099
00100 double* old_a = a; a = new double[nnz];
00101 Array::Copy(old_a,a,nnz);
00102
00103 bLocalCopy = true;
00104 }
00105 }
00106
00107
00108 BOOL SparseMatrixF::IsSamePattern(SparseMatrixF* sm)
00109 {
00110 BOOL isSame = sm->neq == neq
00111 && sm->m_bIsSymmertric == m_bIsSymmertric
00112 && sm->m_eSparseOrientation == m_eSparseOrientation;
00113
00114 if (!isSame)
00115 return FALSE;
00116
00117 ULONG nnz = Nonzeros();
00118 for (ULONG j=0; j<nnz; j++)
00119 if (ci[j] != sm->ci[j])
00120 return false;
00121
00122 for (ULONG j=0;j<neq;j++)
00123 if (adr[j] != sm->adr[j])
00124 return false;
00125
00126 return TRUE;
00127 }
00128
00129 void SparseMatrixF::AddNumbers (double alfa,double* a)
00130 {
00131 ULONG nnz = Nonzeros();
00132 for (ULONG j=0; j<nnz; j++)
00133 this->a[j] += alfa*a[j];
00134 }
00135
00136 void SparseMatrixF::ScaleMatrix(double alfa)
00137 {
00138 ULONG nnz = Nonzeros();
00139 for (ULONG j=0; j<nnz; j++)
00140 this->a[j] *= alfa;
00141 }
00142
00143 long SparseMatrixF::Nonzeros()
00144 {
00145 if (adr==0)
00146 return 0;
00147 if (bJardaConvention)
00148 return (long)(adr[neq]-adr[0]);
00149 else
00150 return (long)(adr[neq-1]);
00151 }
00152
00153 void SparseMatrixF::GetA12block(double* pA12,long c)
00154 {
00155 for (ULONG j=neq-c;j<neq;j++)
00156 {
00157 for (ULONG ad = Adr(j); ad<Adr(j+1); ad++)
00158 {
00159 ULONG i = Ci(ad);
00160 if (i<(neq-c))
00161 pA12[i*c+(j-(neq-c))] = a[ad];
00162 }
00163 }
00164 }
00165
00166 void SparseMatrixF::MulMatrixByVector(double *b,double *c)
00167 {
00168 if (m_bIsSymmertric)
00169 MulSymMatrixByVector(b,c);
00170 else
00171 MulNonsymMatrixByVector(b,c);
00172 }
00173
00174 void SparseMatrixF::MulNonsymMatrixByVector(double *b,double *c)
00175 {
00176 switch (m_eSparseOrientation)
00177 {
00178 case eCompressedColumns:
00179 {
00180 for (ULONG j=0;j<neq;j++)
00181 for (ULONG ad = Adr(j); ad<Adr(j+1); ad++)
00182 {
00183 double A = a[ad];
00184 ULONG i = Ci(ad);
00185 c[i] += A*b[j];
00186 }
00187 }
00188 break;
00189 case eCompressedRows:
00190 {
00191 for (ULONG j=0;j<neq;j++)
00192 for (ULONG ad = Adr(j); ad<Adr(j+1); ad++)
00193 {
00194 double A = a[ad];
00195 ULONG i = Ci(ad);
00196 c[j] += A*b[i];
00197 }
00198 }
00199
00200 break;
00201
00202 default:
00203 break;
00204 }
00205 }
00206
00207 void SparseMatrixF::MulSymMatrixByVector(double *b,double *c)
00208 {
00209 for (ULONG j=0;j<neq;j++)
00210 for (ULONG ad = Adr(j); ad<Adr(j+1); ad++)
00211 {
00212 double A = a[ad];
00213 ULONG i = Ci(ad);
00214 c[i] += A*b[j];
00215 if (i!=j)
00216 c[j] += A*b[i];
00217 }
00218 }
00219
00220 void SparseMatrixF::mxv_scr (double *b,double *c)
00221 {
00222 ULONG i,j,ii,lj,uj;
00223 double s,d;
00224
00225 for (i=0;i<neq;i++){
00226 lj=Adr(i); uj=Adr(i+1);
00227 s=0.0; d=b[i];
00228 for (j=lj;j<uj;j++){
00229 ii=Ci(j);
00230 s+=a[j]*b[ii];
00231 c[ii]+=a[j]*d;
00232 }
00233 c[i]=s;
00234 }
00235 }
00236
00237 void SparseMatrixF::ReadDiagonal(double* dv)
00238 {
00239 for (ULONG j=0;j<neq;j++)
00240 for (ULONG ad = Adr(j); ad<Adr(j+1); ad++)
00241 {
00242 ULONG i = Ci(ad);
00243 if (i==j)
00244 dv[j] = a[ad];
00245 }
00246 }
00247
00248 void SparseMatrixF::LoadMatrix(FILE* stream)
00249 {
00250 neq = 0;adr = NULL;ci = NULL;a = NULL;
00251 bool low = false,up = false;
00252
00253
00254 {
00255
00256 fread(&neq,sizeof(neq),1,stream);
00257
00258 adr = new ULONG[neq+1];
00259 fread(adr,sizeof(ULONG),neq+1,stream);
00260
00261
00262
00263 ci = new ULONG[adr[neq]];
00264 a = new double[adr[neq]];
00265
00266 for (ULONG ad = 0; ad<adr[neq]; ad++)
00267 {
00268
00269 fread(ci+ad,sizeof(ULONG),1,stream);
00270
00271 fread(a+ad,sizeof(double),1,stream);
00272 if (ci[ad]<ad) low = true;
00273 if (ci[ad]>ad) up = true;
00274 }
00275 if (adr[0]==1)
00276 {
00277 Aoffset = 1;
00278 Coffset = 1;
00279 }
00280 bLocalCopy = true;
00281 m_bIsSymmertric = (!low || !up);
00282 }
00283
00284
00285
00286
00287
00288 }
00289
00290 void SparseMatrixF::SaveMatrix(FILE *stream)
00291 {
00292
00293
00294 {
00295
00296
00297 fwrite(&neq,sizeof(neq),1,stream);
00298
00299
00300
00301 fwrite(adr,sizeof(ULONG),neq+1,stream);
00302
00303 for (ULONG ad = 0; ad<adr[neq]; ad++)
00304 {
00305
00306 fwrite(ci+ad,sizeof(ULONG),1,stream);
00307
00308 fwrite(a+ad,sizeof(double),1,stream);
00309 }
00310 }
00311
00312
00313
00314
00315 }
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 DSS_NAMESPASE_END