00001
00002
00003 #include "DSSolver.h"
00004 #include "SparseGridMtxLU.h"
00005 #include "SparseGridMtxLL.h"
00006 #include "SkyLineMtxLDL.h"
00007 #include "MathTracer.h"
00008
00009 DSS_NAMESPASE_BEGIN
00010
00011
00012
00013 DSSolver::DSSolver(MathTracer* pMT)
00014 {
00015 SetMT(pMT);
00016 matrix = NULL;
00017 matrixPD = NULL;
00018 blockSize = 3;
00019 n_blocks = 0;
00020 neq = 0;
00021 m_eState = ISolver::None;
00022 tmpR = NULL;
00023 dom_order = NULL;
00024 mcn = NULL;
00025 MatrixType = eDSSparseMatrix;
00026
00027 orig_matrix = NULL;
00028
00029 fixed = NULL;
00030 lncn = NULL;
00031
00032 SolverType = eDSSFactorizationLDLT;
00033 run_code = 0;
00034 m_bIsSchur = FALSE;
00035
00036 OrderingType = Ordering::ApproxMinimumDegree;
00037
00038 eMT->Writeln("");
00039 eMT->Writeln("___FemCAD--DirectSparse_Solver______");
00040 eMT->Writeln(" Richard Vondracek (c) 2001-2003 ");
00041 }
00042
00043 DSSolver::~DSSolver()
00044 {
00045 Dispose();
00046 }
00047
00048 long DSSolver::Initialize( unsigned char run_code ,eDSSolverType solverType, eDSMatrixType matrixType)
00049 {
00050 Dispose();
00051 eMT->Writeln("");
00052 eMT->Writeln("Sparse Direct initialized!");
00053 this->SolverType = solverType;
00054 this->MatrixType = matrixType;
00055 this->run_code = run_code;
00056
00057 switch (MatrixType)
00058 {
00059 case eDSSparseMatrix:
00060 {
00061 switch(SolverType)
00062 {
00063 case eDSSFactorizationLLTIncomplete:
00064 OrderingType = Ordering::ApproxMinimumDegreeIncomplete;
00065
00066 break;
00067 case eDSSFactorizationLDLTIncomplete:
00068 OrderingType = Ordering::ApproxMinimumDegreeIncomplete;
00069
00070 break;
00071 case eDSSFactorizationLDLT:
00072 case eDSSFactorizationLLT:
00073 case eDSSFactorizationLU:
00074 OrderingType = Ordering::ApproxMinimumDegree;
00075
00076 break;
00077 case eDSSFastCG:
00078 OrderingType = Ordering::None;
00079 break;
00080 default:
00081 OrderingType = Ordering::None;
00082 break;
00083 }
00084 break;
00085 }
00086 case eDSSkylineMatrix:
00087 {
00088 switch(SolverType)
00089 {
00090 case eDSSFactorizationLLTIncomplete:
00091 OrderingType = Ordering::ReverseCuthillMcKee;
00092 break;
00093 default:
00094 OrderingType = Ordering::None;
00095 break;
00096 }
00097 }
00098 break;
00099 }
00100
00101 m_eState = ISolver::Initialized;
00102 return 1;
00103 }
00104
00105 void DSSolver::SetMT(MathTracer* pMT)
00106 {
00107 if (pMT) eMT = pMT; else eMT = &MT;
00108 }
00109
00110 void DSSolver::Dispose()
00111 {
00112 if (matrix) { delete matrix;matrix = NULL; }
00113 if (matrixPD) { delete matrixPD;matrixPD = NULL; }
00114 if (orig_matrix){ delete orig_matrix;orig_matrix = NULL;}
00115 if (tmpR ) { delete [] tmpR;tmpR = NULL; }
00116 if (dom_order) { delete dom_order;dom_order = NULL; }
00117 if (mcn) { delete mcn;mcn = NULL; }
00118
00119 if (fixed) { delete fixed;fixed = NULL;}
00120 if (lncn) { delete lncn;lncn = NULL;}
00121 }
00122
00123 BOOL DSSolver::SetOrderingType(Ordering::Type otype)
00124 {
00125 switch(otype)
00126 {
00127 case Ordering::None:
00128 case Ordering::MinimumDegree:
00129 case Ordering::ApproxMinimumDegree:
00130 case Ordering::ReverseCuthillMcKee:
00131 case Ordering::CuthillMcKee:
00132 case Ordering::NestedGraphBisection:
00133 #ifdef _LINK_METIS_
00134 case Ordering::MetisND:
00135 #endif
00136 #ifdef _LINK_COLAMD_
00137 case Ordering::ColAMD:
00138 #endif
00139 OrderingType = otype;
00140 return TRUE;
00141 default:
00142 eMT->Writeln("Selected ordering type is not implemented! Default value is used.");
00143 break;
00144 }
00145 return FALSE;
00146 }
00147
00148 BOOL DSSolver::LoadMatrix(unsigned long neq,unsigned char block_size,double * a,unsigned long * ci,unsigned long * adr )
00149 {
00150 SparseMatrixF loc_sm(neq,a,ci,adr);
00151 return LoadMatrix(&loc_sm,block_size);
00152 }
00153
00154 BOOL DSSolver::LoadMatrix(SparseMatrixF* smt,unsigned char block_size)
00155 {
00156 SetMatrixPattern(smt,block_size);
00157 m_sm.CreateLocalCopy();
00158 return 1;
00159 }
00160
00161 BOOL DSSolver::SetMatrixPattern(SparseMatrixF* smt,unsigned char block_size)
00162 {
00163 Dispose();
00164 m_sm = *smt;
00165 blockSize = block_size;
00166 neq = (long)m_sm.neq;
00167
00168 eMT->Writeln("Loading sparse matrix");
00169 sprintf(str," number of equations : %ld ",m_sm.neq);
00170 eMT->Writeln(str);
00171 sprintf(str," number of nonzeros : %ld nonzeros",m_sm.Nonzeros());
00172 eMT->Writeln(str);
00173
00174 if (SolverType==eDSSDiagonalScaling)
00175 {
00176 matrixD.Init(neq);
00177 m_sm.ReadDiagonal(matrixD.DataPtr());
00178 }
00179 m_eState = ISolver::Initialized;
00180 return 1;
00181 }
00182
00183 BOOL DSSolver::IsInitialized()
00184 {
00185 return m_eState != ISolver::None;
00186 }
00187
00188 BOOL DSSolver::IsFactorized()
00189 {
00190 return m_eState == ISolver::Factorized;
00191 }
00192
00193 BOOL DSSolver::IsAllocated()
00194 {
00195 return m_eState == ISolver::Allocated
00196 || m_eState == ISolver::Factorized;
00197 }
00198
00199 BOOL DSSolver::IsSchur()
00200 {
00201 return m_bIsSchur;
00202 }
00203
00204 SparseGridMtx* DSSolver::CreateNewSparseGridMtx(IntArrayList* fixed)
00205 {
00206 long no_nonzeros = m_sm.Nonzeros();
00207 eMT->CS();
00208
00209 SparseConectivityMtxII* block_conmtx = new SparseConectivityMtxII(m_sm,this->mcn,blockSize);
00210 block_conmtx->eMT = eMT;
00211
00212 long no_init_blocks = (block_conmtx->Nonzeros())/2+block_conmtx->N();
00213 n_blocks = block_conmtx->N();
00214
00215 Ordering* order = block_conmtx->GetPermutationAndPattern(OrderingType,fixed);
00216 delete block_conmtx;
00217 eMT->Writeln(eMT->MC_());
00218 block_conmtx = NULL;
00219
00220 SparseGridMtx* matrix = NULL;
00221
00222 eMT->Writeln("Allocating block sparse matrix");
00223 switch(SolverType)
00224 {
00225 case eDSSFactorizationLDLTIncomplete:
00226 matrix = new SparseGridMtxLDL(m_sm,blockSize,order,mcn,eMT,TRUE);
00227 break;
00228 case eDSSFactorizationLLTIncomplete:
00229 matrix = new SparseGridMtxLL(m_sm,blockSize,order,mcn,eMT,TRUE);
00230 break;
00231 case eDSSFactorizationLDLT:
00232 matrix = new SparseGridMtxLDL(m_sm,blockSize,order,mcn,eMT,TRUE);
00233 break;
00234 case eDSSFactorizationLLT:
00235 matrix = new SparseGridMtxLL(m_sm,blockSize,order,mcn,eMT,TRUE);
00236 break;
00237 case eDSSFactorizationLU:
00238 matrix = new SparseGridMtxLU(m_sm,blockSize,order,mcn,eMT,TRUE);
00239 break;
00240 case eDSSFastCG:
00241 matrix = new SparseGridMtxLDL(m_sm,blockSize,order,mcn,eMT,TRUE);
00242 break;
00243 default:
00244 eMT->Writeln("Unknown solver type.");
00245 break;
00246 }
00247
00248 if (matrix==NULL)
00249 return NULL;
00250
00251 order = NULL;
00252
00253 matrix->WriteStatistics(no_init_blocks,no_nonzeros);
00254
00255 if (fixed!=NULL)
00256 {
00257 sprintf(str," Schur complement : %ldx%ld\n",fixed->Count,fixed->Count);
00258 eMT->Write(str);
00259 }
00260
00261 if (tmpR) {delete [] tmpR;tmpR = NULL;}
00262 tmpR = new double[n_blocks*blockSize];
00263 Array::Clear(tmpR,0,n_blocks*blockSize);
00264 return matrix;
00265 }
00266
00267 BOOL DSSolver::PreFactorize()
00268 {
00269 if (m_sm.neq == 0 || neq == 0)
00270 {
00271 eMT->Write("Can't factorize empty matrix!");
00272 return FALSE;
00273 }
00274
00275
00276 delete matrix;
00277 matrix = NULL;
00278 delete matrixPD;
00279 matrixPD = NULL;
00280
00281 switch (MatrixType)
00282 {
00283 case eDSSparseMatrix:
00284 {
00285 if (IsSchur())
00286 {
00287 return PreFactorizeSchur();
00288 }
00289 else
00290 {
00291 this->matrix = CreateNewSparseGridMtx();
00292 if (this->matrix==NULL)
00293 return FALSE;
00294 }
00295 }
00296 break;
00297
00298 case eDSSkylineMatrix:
00299 switch (SolverType)
00300 {
00301 case eDSSFactorizationLDLT:
00302
00303
00304 case eDSSFactorizationLDLTIncomplete:
00305 case eDSSFactorizationLLTIncomplete:
00306 case eDSSFactorizationLLT:
00307 case eDSSFactorizationLU:
00308 case eDSSFastCG:
00309 eMT->Writeln("Skyline matrix doesn't support this solver type!");
00310 break;
00311 default:
00312 eMT->Writeln("Unknown solver type.");
00313 break;
00314 }
00315 break;
00316 default:
00317 eMT->Writeln("Unknown matrix storage type.");
00318 break;
00319 }
00320
00321 if (this->matrix==NULL)
00322 return FALSE;
00323
00324 m_eState = ISolver::Allocated;
00325 return TRUE;
00326 }
00327
00328 void DSSolver::LoadZeros()
00329 {
00330 if (matrix==NULL){
00331 eMT->Writeln("First allocate the matrix!");
00332 return;
00333 }
00334 if (matrix)
00335 matrix->LoadZeros();
00336 m_eState = ISolver::Allocated;
00337 }
00338
00339 double& DSSolver::ElementAt(int i, int j)
00340 {
00341 return matrix->ElementAt(i,j);
00342 }
00343
00344 BOOL DSSolver::LoadNumbers (SparseMatrixF* sm)
00345 {
00346 if (matrix==NULL){
00347 eMT->Writeln("First load matrix.");
00348 return 0;
00349 }
00350 if (!m_sm.IsSamePattern(sm))
00351 {
00352 eMT->Writeln("Matrix has different nonzero pattern!");
00353 return 0;
00354 }
00355 if (m_sm.a==NULL || m_sm.a!=sm->a)
00356 {
00357 m_sm = *sm;
00358 }
00359
00360 eMT->CS();eMT->Write("Loading matrix data : ");
00361 matrix->LoadMatrixNumbers(*sm);
00362 eMT->Writeln(eMT->MC_());
00363
00364 m_eState = ISolver::Allocated;
00365 return 1;
00366 }
00367
00368 BOOL DSSolver::AddNumbers (double alfa,SparseMatrixF* smtx)
00369 {
00370 if (m_eState != ISolver::Allocated)
00371 {
00372 eMT->Writeln("First load matrix.");
00373 return 0;
00374 }
00375
00376 m_sm.AddNumbers(alfa,smtx->a);
00377 return TRUE;
00378 }
00379
00380 BOOL DSSolver::ScaleMatrix(double alfa)
00381 {
00382 if (m_eState != ISolver::Allocated)
00383 {
00384 eMT->Writeln("First load matrix.");
00385 return 0;
00386 }
00387
00388 m_sm.ScaleMatrix(alfa);
00389 return TRUE;
00390 }
00391
00392 SparseMatrixF* DSSolver::GetSparseMatrix()
00393 {
00394 if (m_eState != ISolver::Allocated)
00395 {
00396 eMT->Writeln("First load matrix.");
00397 return NULL;
00398 }
00399 return &m_sm;
00400 }
00401
00402
00403 void DSSolver::WriteFactorizationInfo()
00404 {
00405 switch(SolverType)
00406 {
00407 case eDSSFactorizationLDLTIncomplete:
00408 eMT->Write("Incomplete LDL factorization : ");
00409 break;
00410 case eDSSFactorizationLLTIncomplete:
00411 eMT->Write("Incomplete LL factorization : ");
00412 break;
00413 case eDSSFactorizationLDLT:
00414 eMT->Write("Numerical LDL factorization : ");
00415 break;
00416 case eDSSFactorizationLLT:
00417 eMT->Write("Numerical LLT factorization : ");
00418 break;
00419 case eDSSFactorizationLU:
00420 eMT->Write("Numerical LU factorization : ");
00421 break;
00422 case eDSSFastCG:
00423 break;
00424 default:
00425 eMT->Writeln("Unknown solver type.");
00426 break;
00427 }
00428 }
00429
00430 BOOL DSSolver::ReFactorize( )
00431 {
00432 if (matrix==NULL && matrixPD==NULL){
00433 eMT->Writeln("The matrix has to be loaded prior to the factorization.");
00434 return FALSE;
00435 }
00436
00437 if (!IsAllocated()){
00438 eMT->Writeln("The matrix has to be allocated (PreFactorize) prior to the factorization.");
00439 return FALSE;
00440 }
00441 eMT->CS();
00442
00443 DenseMatrixArithmetics::zero_pivots = 0;
00444
00445 WriteFactorizationInfo();
00446 switch(SolverType)
00447 {
00448 case eDSSFactorizationLDLTIncomplete:
00449 case eDSSFactorizationLLTIncomplete:
00450 case eDSSFactorizationLDLT:
00451 case eDSSFactorizationLLT:
00452 case eDSSFactorizationLU:
00453 if (matrixPD!=NULL)
00454 matrixPD->SchurComplementFactorization();
00455 else
00456 matrix->Factorize();
00457 break;
00458 case eDSSFastCG:
00459 break;
00460 default:
00461 eMT->Writeln("Unknown solver type.");
00462 break;
00463 }
00464
00465 eMT->Writeln(eMT->MC_());
00466 if (DenseMatrixArithmetics::zero_pivots>0)
00467 {
00468 sprintf(str,"Warning: %ld zero pivot",DenseMatrixArithmetics::zero_pivots);
00469 eMT->Write(str);
00470 if (DenseMatrixArithmetics::zero_pivots == 1)
00471 eMT->Writeln("has been found during the factorization !!!");
00472 else
00473 eMT->Writeln("s have been found during the factorization !!!");
00474 }
00475
00476 m_eState = ISolver::Factorized;
00477 return TRUE;
00478 }
00479
00480
00481 BOOL DSSolver::Solve(double * r, double * f )
00482 {
00483 if (!IsFactorized())
00484 {
00485 if (!Factorize())
00486 return FALSE;
00487 }
00488
00489 if (matrix==NULL)
00490 {
00491 if (matrixD.N())
00492 matrixD.DiagonalSolve(f,r);
00493 else
00494 if (r!=f)
00495 Array::Copy(f,r,neq);
00496 return FALSE;
00497 }
00498
00499 if (mcn)
00500 {
00501 long i;
00502 long* perm = mcn->perm->Items;
00503 for (i=0; i<neq; i++)
00504 tmpR[perm[i]] = f[i];
00505
00506 matrix->Solve(tmpR,tmpR);
00507
00508 for (i=0; i<neq; i++)
00509 r[i] = tmpR[perm[i]];
00510 }
00511 else
00512 {
00513 if (neq%blockSize)
00514 {
00515 Array::Copy(f,tmpR,neq);
00516 matrix->Solve(tmpR,tmpR);
00517 Array::Copy(tmpR,r,neq);
00518 }
00519 else
00520 {
00521 matrix->Solve(f,r);
00522 }
00523 }
00524 return TRUE;
00525 }
00526
00527
00528 long DSSolver::Close ( )
00529 {
00530 Dispose();
00531 return 0;
00532 }
00533
00534 void DSSolver::StartSolverWriteInfo()
00535 {
00536 switch (SolverType)
00537 {
00538 case eDSSFactorizationLDLT:
00539 eMT->Writeln("Sparse matrix LDL factorization");
00540 break;
00541 case eDSSFactorizationLLT:
00542 eMT->Writeln("Sparse matrix LL factorization");
00543 break;
00544 case eDSSFactorizationLU:
00545 eMT->Writeln("Sparse matrix LU factorization");
00546 break;
00547 case eDSSFactorizationLDLTIncomplete:
00548 eMT->Writeln("Incomplete matrix LDL factorization");
00549 break;
00550 case eDSSFactorizationLLTIncomplete:
00551 eMT->Writeln("Incomplete matrix LL factorization");
00552 break;
00553 default:
00554 break;
00555 }
00556 sprintf(str,"Starting on : %s",eMT->NowString());
00557 eMT->Write(str);
00558 }
00559
00560 void DSSolver::EndSolverWriteInfo()
00561 {
00562 eMT->Write( "Ending on : ");
00563 eMT->Write(eMT->NowString());
00564 if (matrix){
00565 sprintf(str,"Computational effort : %ld block multiplications",matrix->No_Multiplications());
00566 }
00567 eMT->Writeln(str);
00568 }
00569
00570 BOOL DSSolver::Factorize()
00571 {
00572 if (m_sm.neq == 0)
00573 {
00574 eMT->Write("Can't factorize empty matrix!");
00575 return FALSE;
00576 }
00577
00578 long usedmem = 0;
00579 clock_t solution_start = clock();
00580
00581 StartSolverWriteInfo();
00582 if (!IsAllocated())
00583 {
00584 if (!PreFactorize( ))
00585 return FALSE;
00586 }
00587 if (!ReFactorize())
00588 return FALSE;
00589
00590 clock_t solution_end = clock();
00591 EndSolverWriteInfo();
00592
00593 double solution_duration = (double)(solution_end-solution_start) / CLOCKS_PER_SEC;
00594 sprintf(str,"Whole solution (end-start) : %0.3f s",solution_duration);
00595 eMT->Writeln(str);
00596
00597 usedmem = 0;
00598
00599
00600 return TRUE;
00601 }
00602
00603 void DSSolver::SetSM(SparseMatrixF* sm)
00604 {
00605 this->m_sm = *sm;
00606 }
00607
00608 double DSSolver::GetFactorizationError()
00609 {
00610 if (neq==0)
00611 return 0.0;
00612 if (m_sm.a==NULL)
00613 return -1;
00614
00615 LargeVector r(neq);
00616 LargeVector b(neq);
00617 LargeVector x(neq);
00618 for (int i=0; i<neq; i++)
00619 b.DataPtr()[i] = 1.0;
00620
00621 Solve(x.DataPtr(),b.DataPtr());
00622
00623
00624 MulMatrixByVector(x.DataPtr(),r.DataPtr());
00625 r.LinComb(-1,r,b);
00626
00627
00628 double err = LargeVector::InnerProduct(r.DataPtr(),r.DataPtr(),neq);
00629 return sqrt(err)/neq;
00630 }
00631
00632
00633
00634
00635 void DSSolver::ExpandMCN(IntArrayList& mcn)
00636 {
00637 long nb = mcn.Count/blockSize;
00638 for (long b=0; b<nb; b++)
00639 {
00640 bool fixed = false;
00641 bool nonfixed = false;
00642 for (long i=0; i<blockSize; i++)
00643 {
00644 long v = mcn[b*blockSize+i];
00645 if (v<0) fixed=true;
00646 if (v>=0) nonfixed=true;
00647 }
00648 if (fixed && nonfixed)
00649 {
00650 for (long i=0; i<blockSize; i++)
00651 {
00652 long cfix = -1;
00653 long v = mcn[b*blockSize+i];
00654 if (v<0)
00655 {
00656 cfix=v;
00657 mcn[b*blockSize+i] = -1;
00658 }
00659 mcn.Add(cfix);
00660 }
00661 }
00662 }
00663 }
00664
00665
00666
00667
00668
00669 BOOL DSSolver::LoadMCN(ULONG n_blocks,unsigned char block_size,long * mcn,BOOL bIsSchur)
00670 {
00671 this->blockSize = block_size;
00672 IntArrayList* mcn_order = new IntArrayList(n_blocks*block_size);
00673 mcn_order->Alloc();
00674 Array::Copy(mcn,mcn_order->Items,n_blocks*block_size);
00675 m_bIsSchur = bIsSchur;
00676 return LoadMCN_int(mcn_order);
00677 }
00678
00679 BOOL DSSolver::LoadMCN_int(IntArrayList* mcn_order)
00680 {
00681 ExpandMCN(*mcn_order);
00682 this->n_blocks = mcn_order->Count/blockSize;
00683
00684 long i=0;
00685 IntArrayList* mcn_perm = new IntArrayList(mcn_order->Count);
00686 mcn_perm->Alloc();
00687
00688 long no_noncondensed_DOFs = 0;
00689
00690 for (i=0; i<mcn_perm->Count; i++)
00691 mcn_perm->Items[i] = -1;
00692
00693 for (i=0; i<mcn_perm->Count; i++)
00694 {
00695 long oi = mcn_order->Items[i];
00696 if (oi>=0)
00697 mcn_perm->Items[oi] = i;
00698 else if (oi==-1)
00699 {
00700
00701 } else {
00702 mcn_perm->Items[-(oi+2)] = i;
00703 no_noncondensed_DOFs++;
00704 }
00705 }
00706
00707 if (this->mcn) delete this->mcn;
00708 this->mcn = new Ordering(mcn_perm,mcn_order);
00709
00710 return CreateFixedArray(no_noncondensed_DOFs);
00711 }
00712
00713 BOOL DSSolver::LoadMCN (IntArrayList& mcn)
00714 {
00715 IntArrayList* mcn_new = new IntArrayList(((mcn.Count-1)/blockSize+1)*blockSize);
00716 mcn_new->Alloc();
00717 mcn_new->Initialize(-1);
00718 Array::Copy(mcn.Items,0,mcn_new->Items,0,mcn.Count);
00719 return LoadMCN_int(mcn_new);
00720 }
00721
00722 BOOL DSSolver::CreateFixedArray(long no_noncondensed_DOFs)
00723 {
00724 fixed = new IntArrayList();
00725 lncn = new IntArrayList(no_noncondensed_DOFs);
00726
00727 lncn->Alloc();
00728 long lastFixed = -1;
00729
00730 long i=0;
00731 for (i=0; i<mcn->perm->Count; i++)
00732 {
00733 long oi = mcn->order->Items[i];
00734 if (oi<-1)
00735 {
00736 fixed->AddIfNotLast(i/blockSize);
00737
00738
00739
00740 long cn = -(oi+2);
00741 if (cn<lastFixed) {
00742 eMT->Writeln("Error: Array of fixed nodes is not monotonously growing !!");
00743 m_eState = ErrorInMCN;
00744 } else
00745 if (cn>=neq) {
00746 eMT->Writeln("Error: Array of fixed nodes contains index higher than (neq) !!");
00747 m_eState = ErrorInMCN;
00748 } else if (cn<(neq-no_noncondensed_DOFs)) {
00749 eMT->Writeln("Error: Array of fixed nodes contains index lower than (neq-no_fixed_DOFs) !!");
00750 m_eState = ErrorInMCN;
00751 } else
00752 lncn->Items[cn-(neq-no_noncondensed_DOFs)] = lastFixed = cn;
00753 }
00754 }
00755
00756
00757
00758
00759
00760
00761
00762 return m_eState != ErrorInMCN;
00763 }
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839 BOOL DSSolver::PreFactorizeSchur()
00840 {
00841 if (mcn == 0)
00842 {
00843 eMT->Writeln("The MCN vector was not set! It is needed for Schur complement condensation.");
00844 return FALSE;
00845 }
00846
00847 StoreFixedLastPermutation_dom_order();
00848
00849 delete matrixPD;
00850 matrixPD = NULL;
00851
00852 SparseGridMtx* matrix = CreateNewSparseGridMtx(fixed);
00853 matrixPD = new SparseGridMtxPD(matrix,fixed->Count);
00854 m_eState = ISolver::Allocated;
00855
00856 return TRUE;
00857 }
00858
00859
00860
00861 void DSSolver::StoreFixedLastPermutation_dom_order()
00862 {
00863
00864 if (dom_order) delete dom_order;
00865 dom_order = new IntArrayList(m_sm.neq);
00866 dom_order->Alloc();
00867
00868 int lastdo = 0;
00869
00870 long* dmap = new long[neq];
00871 Array::Clear(dmap,0,neq);
00872
00873 if (lncn->Count)
00874 {
00875 long*pI = lncn->Items;
00876 long idx = lncn->Count-1;
00877 for(long* p = pI+idx; p>=pI; p--,idx--)
00878 dmap[*p] = ~(idx+(neq-lncn->Count));
00879 }
00880
00881 for (long d=0; d<neq; d++)
00882 {
00883 if (dmap[d]==0)
00884 {
00885 dmap[d] = ~lastdo;
00886 dom_order->Items[lastdo++] = d;
00887 }
00888 }
00889 Array::Copy(lncn->Items,dom_order->Items+neq-lncn->Count,lncn->Count);
00890 delete dmap; dmap = NULL;
00891 }
00892
00893 void DSSolver::condense(double *a,double *lhs,double *rhs,long tc)
00894 {
00895 if (!IsSchur())
00896 {
00897 if (m_eState == ISolver::Initialized)
00898 {
00899 m_bIsSchur = TRUE;
00900 }
00901 else
00902 {
00903 eMT->Writeln("Error : If you want to call condense you should switch IsSchur to TRUE in LoadMCN!");
00904 return;
00905 }
00906 }
00907
00908 if (tc == 1)
00909 {
00910 if (m_eState!=ISolver::Factorized)
00911 if (!Factorize())
00912 return;
00913
00914 matrixPD->WriteCondensedMatrixA22(a,mcn,lncn);
00915
00916 long i=0;
00917 long* perm = mcn->perm->Items;
00918 long* dor = dom_order->Items;
00919 for (i=0; i<neq; i++)
00920 tmpR[perm[dor[i]]] = rhs[i];
00921
00922 matrixPD->Sub_A21_A11inv(tmpR);
00923
00924 for (i=0; i<neq; i++)
00925 rhs[i] = tmpR[perm[dor[i]]];
00926 }
00927 else if (tc==2)
00928 {
00929 long i;
00930 long* perm = mcn->perm->Items;
00931 long* dor = dom_order->Items;
00932 for (i=0; i<neq-lncn->Count; i++)
00933 tmpR[perm[dor[i]]] = rhs[i];
00934
00935 for (i=neq-lncn->Count; i<neq; i++)
00936 tmpR[perm[dor[i]]] = lhs[i];
00937
00938 matrixPD->Sub_A11inv_A12(tmpR);
00939
00940 for (i=0; i<neq; i++)
00941 lhs[i] = tmpR[perm[dor[i]]];
00942 }
00943 else if (tc==3)
00944 {
00945 long i;
00946 long* perm = mcn->perm->Items;
00947 long* dor = dom_order->Items;
00948 for (i=0; i<neq-lncn->Count; i++)
00949 tmpR[perm[dor[i]]] = rhs[i];
00950
00951 matrixPD->SolveA11(tmpR);
00952
00953 for (i=0; i<neq-lncn->Count; i++)
00954 lhs[i] = tmpR[perm[dor[i]]];
00955 }
00956 }
00957
00958 void DSSolver::GetA12block(double *pA12)
00959 {
00960 m_sm.GetA12block(pA12,lncn->Count);
00961 }
00962
00963
00964
00965 void DSSolver::MulMatrixByVector(double *b, double *c)
00966 {
00967 if (SolverType==eDSSFastCG)
00968 {
00969 LargeVectorAttach vb(neq,b);
00970 LargeVectorAttach vx(neq,c);
00971 matrix->MultiplyByVector(vb,vx);
00972 }
00973
00974 if (SolverType==eDSSFactorizationLU)
00975 m_sm.MulNonsymMatrixByVector(b,c);
00976 else
00977 m_sm.MulMatrixByVector(b,c);
00978
00979 }
00980
00981
00982
00983
00984 int DSSolver::PreCG(double* b_, double* x_,double epsilon,int max_iter)
00985 {
00986 double epsilon2 = epsilon*epsilon;
00987
00988 LargeVectorAttach b(neq,b_);
00989 LargeVectorAttach x(neq,x_);
00990
00991
00992 LargeVector r(neq);
00993 MulMatrixByVector(x.DataPtr(),r.DataPtr());
00994
00995
00996 r.LinComb(-1,r,b);
00997
00998
00999 LargeVector d(neq);
01000 Solve(d.DataPtr(),r.DataPtr());
01001
01002
01003 double delta_new = LargeVector::InnerProduct(r.DataPtr(),d.DataPtr(),neq);
01004
01005 if (delta_new==0)
01006 {
01007 eMT->Writeln("Zero right hand side vector.");
01008 return 0;
01009 }
01010 double delta_0 = delta_new;
01011
01012 eMT->Writeln("Starting PCG ");
01013 eMT->Write(" epsilon : ");eMT->Write(epsilon);eMT->Writeln();
01014 eMT->Write(" residual : ");eMT->Writeln();
01015
01016 LargeVector q(neq);
01017 LargeVector s(neq);
01018 double alfa = 0;
01019 int iter;
01020 for (iter=1;iter<max_iter; iter++)
01021 {
01022
01023 MulMatrixByVector(d.DataPtr(),q.DataPtr());
01024
01025
01026 double a1 = LargeVector::InnerProduct(d.DataPtr(),q.DataPtr(),neq);
01027 alfa = delta_new/a1;
01028
01029
01030 x.AddMult(alfa,d);
01031
01032
01033 if (iter%50 == 0)
01034 {
01035
01036 MulMatrixByVector(x.DataPtr(),r.DataPtr());
01037 r.LinComb(-1,r,b);
01038
01039
01040 }
01041 else
01042 {
01043
01044 r.AddMult(-alfa,q);
01045 }
01046
01047
01048 Solve(s.DataPtr(),r.DataPtr());
01049
01050 double delta_old = delta_new;
01051
01052 delta_new = LargeVector::InnerProduct(r.DataPtr(),s.DataPtr(),neq);
01053
01054 double beta = delta_new/delta_old;
01055
01056
01057
01058
01059 d.LinComb(beta,d,s);
01060
01061 eMT->Write(" ");
01062 eMT->Write(sqrt(delta_new/delta_0));
01063 eMT->Writeln();
01064
01065
01066 if (fabs(delta_new) < fabs(epsilon2*delta_0))
01067 {
01068 eMT->Write("Success after ");
01069 eMT->Write(iter);
01070 eMT->Writeln(" iterations");
01071 return iter;
01072 }
01073 }
01074 eMT->Write("Failed after ");
01075 eMT->Write(iter);
01076 eMT->Writeln(" iterations");
01077 return -iter;
01078 }
01079
01080
01081
01082 int DSSolver::CG(double* b_, double* x_,double epsilon,int max_iter)
01083 {
01084 double epsilon2 = epsilon*epsilon;
01085
01086 LargeVectorAttach b(neq,b_);
01087 LargeVectorAttach x(neq,x_);
01088
01089
01090 LargeVector r(neq);
01091 MulMatrixByVector(x.DataPtr(),r.DataPtr());
01092
01093
01094
01095 r.LinComb(-1,r,b);
01096
01097
01098 LargeVector d(r);
01099
01100
01101 double delta_new = LargeVector::InnerProduct(r.DataPtr(),r.DataPtr(),neq);
01102
01103 if (delta_new==0)
01104 {
01105 eMT->Writeln("Zero right hand side vector.");
01106
01107 }
01108 double delta_0 = delta_new;
01109
01110 eMT->Writeln("Starting PCG ");
01111 eMT->Write(" epsilon : ");eMT->Write(epsilon);eMT->Writeln();
01112 eMT->Write(" residual : ");eMT->Writeln();
01113
01114 LargeVector q(neq);
01115 LargeVector s(neq);
01116 double alfa = 0;
01117 int iter;
01118 for (iter=1;iter<max_iter; iter++)
01119 {
01120
01121 MulMatrixByVector(d.DataPtr(),q.DataPtr());
01122
01123
01124 double a1 = LargeVector::InnerProduct(d.DataPtr(),q.DataPtr(),neq);
01125 alfa = delta_new/a1;
01126
01127
01128 x.AddMult(alfa,d);
01129
01130
01131 if (iter%50 == 0)
01132 {
01133
01134 MulMatrixByVector(x.DataPtr(),r.DataPtr());
01135 r.LinComb(-1,r,b);
01136
01137
01138 }
01139 else
01140 {
01141
01142 r.AddMult(-alfa,q);
01143 }
01144
01145 double delta_old = delta_new;
01146
01147 delta_new = LargeVector::InnerProduct(r.DataPtr(),r.DataPtr(),neq);
01148
01149 double beta = delta_new/delta_old;
01150 double conv = delta_new/delta_0;
01151
01152 eMT->Write(" ");
01153 eMT->Write(sqrt(conv));
01154
01155 eMT->Writeln();
01156
01157
01158
01159
01160 d.LinComb(beta,d,r);
01161
01162
01163 if (fabs(conv) < fabs(epsilon2))
01164 {
01165 return iter;
01166 }
01167 }
01168 return -iter;
01169 }
01170
01171
01172
01173
01174 DSS_NAMESPASE_END