00001
00002
00003 #include "SparseGridMtxLL.h"
00004
00005 DSS_NAMESPASE_BEGIN
00006
00007
00008
00009 SparseGridMtxLL::SparseGridMtxLL(SparseMatrixF& sm,BYTE block_size,Ordering* block_order,MathTracer* eMT,BOOL load_data)
00010 : SparseGridMtx(sm,block_size,block_order,eMT)
00011 {
00012 this->BlockArith->prefered_decomposition = eLL_decomposition;
00013 IConectMatrix* bskl = block_order->cm;
00014
00015 this->AlocateMemoryByPattern(bskl);
00016 ComputeBlocks();
00017
00018 if (sm.a!=NULL && load_data!=FALSE)
00019 LoadMatrixNumbers(sm);
00020 }
00021
00022
00023
00024 SparseGridMtxLL::SparseGridMtxLL(SparseMatrixF& sm,BYTE block_size,Ordering* block_order,Ordering* node_order,MathTracer* eMT,BOOL load_data )
00025 : SparseGridMtx(sm,block_size,block_order,node_order,eMT)
00026 {
00027 this->BlockArith->prefered_decomposition = eLL_decomposition;
00028 IConectMatrix* bskl = block_order->cm;
00029
00030 this->AlocateMemoryByPattern(bskl);
00031 ComputeBlocks();
00032
00033 if (sm.a!=NULL && load_data!=FALSE)
00034 LoadMatrixNumbers(sm);
00035 }
00036
00037 double& SparseGridMtxLL::ElementAt(int i, int j)
00038 {
00039 long* blockP = this->block_order->perm->Items;
00040 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00041
00042 long aux_bi_idx = 0;
00043 long aux_bj_idx = 0;
00044
00045 long nj = (nodeP==NULL)?j:nodeP[j];
00046 long sj = nj % block_size;
00047 long bj = nj / block_size;
00048 long nbj = (blockP==NULL)?bj:blockP[bj];
00049
00050 SparseGridColumn* Column_nbj = Columns[nbj];
00051
00052 double* column_nbj_bi_data = NULL;
00053 double* column_nbi_bj_data = NULL;
00054
00055 long ni = (nodeP==NULL)?i:nodeP[i];
00056 long si = ni % block_size;
00057 long bi = ni / block_size;
00058 long nbi = (blockP==NULL)?bi:blockP[bi];
00059
00060 SparseGridColumn* Column_nbi = Columns[nbi];
00061
00062 aux_bj_idx = aux_bi_idx = -1;
00063 if (nbi<nbj)
00064 {
00065 aux_bi_idx = Column_nbj->FindExistingBlockIndex(nbi);
00066 column_nbj_bi_data = Columns_data+Column_nbj->column_start_idx+aux_bi_idx*block_storage;
00067 }
00068 else if (nbi>nbj)
00069 {
00070 aux_bj_idx = Column_nbi->FindExistingBlockIndex(nbj);
00071 column_nbi_bj_data = Columns_data+Column_nbi->column_start_idx+aux_bj_idx*block_storage;
00072 }
00073
00074 if (aux_bi_idx<0 && aux_bj_idx<0)
00075 {
00076 if (si<sj)
00077 return Columns_data[nbi*block_storage+si+sj*block_size];
00078 else
00079 return Columns_data[nbi*block_storage+sj+si*block_size];
00080 }
00081 else if (aux_bi_idx>=0)
00082 return column_nbj_bi_data[block_size*sj + si];
00083 else
00084 return column_nbi_bj_data[block_size*si + sj];
00085 }
00086
00087 void SparseGridMtxLL::LoadZeros()
00088 {
00089 memset(this->Columns_data,0,columns_data_length*sizeof(double));
00090 for (long bj=0; bj<n_blocks; bj++)
00091 for (long j=0; j<block_size; j++)
00092 Columns_data[bj*block_storage + block_size*j + j] = 1.0;
00093
00094 long* blockP = this->block_order->perm->Items;
00095 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00096
00097 long neq = n_blocks*block_size - noDummyDOFs;
00098 for (long i=0; i<neq; i++)
00099 {
00100 long ni = (nodeP==NULL)?i:nodeP[i];
00101 long si = ni % block_size;
00102 long bi = ni / block_size;
00103 long nbi = (blockP==NULL)?bi:blockP[bi];
00104 Columns_data[nbi*block_storage+si+si*block_size] = 0.0;
00105 }
00106 }
00107
00108 void SparseGridMtxLL::LoadMatrixNumbers(SparseMatrixF& sm)
00109 {
00110 LoadZeros();
00111 long* blockP = this->block_order->perm->Items;
00112 long* nodeP = this->node_order?node_order->perm->Items:NULL;
00113
00114 long aux_bi_idx = 0;
00115 long aux_bj_idx = 0;
00116
00117 for (ULONG j=0; j<sm.neq; j++)
00118 {
00119 long nj = (nodeP==NULL)?j:nodeP[j];
00120
00121 long sj = nj % block_size;
00122 long bj = nj / block_size;
00123 long nbj = (blockP==NULL)?bj:blockP[bj];
00124
00125 SparseGridColumn* Column_nbj = Columns[nbj];
00126 SparseGridColumn* Column_nbi = NULL;
00127
00128 double* column_nbj_bi_data = NULL;
00129 double* column_nbi_bj_data = NULL;
00130
00131 long old_nbi = -1;
00132
00133 for (ULONG ad = sm.Adr(j); ad<sm.Adr(j+1); ad++)
00134 {
00135 long i = (long)sm.Ci(ad);
00136 long ni = (nodeP==NULL)?i:nodeP[i];
00137 double val = sm.a[ad];
00138
00139 long si = ni % block_size;
00140 long bi = ni / block_size;
00141 long nbi = (blockP==NULL)?bi:blockP[bi];
00142
00143 if (old_nbi!=nbi)
00144 {
00145 Column_nbi = Columns[nbi];
00146 aux_bj_idx = aux_bi_idx = -1;
00147 if (nbi<nbj)
00148 {
00149 aux_bi_idx = Column_nbj->FindExistingBlockIndex(nbi);
00150 column_nbj_bi_data = Columns_data+Column_nbj->column_start_idx+aux_bi_idx*block_storage;
00151 }
00152 else if (nbi>nbj)
00153 {
00154 aux_bj_idx = Column_nbi->FindExistingBlockIndex(nbj);
00155 column_nbi_bj_data = Columns_data+Column_nbi->column_start_idx+aux_bj_idx*block_storage;
00156 }
00157 old_nbi=nbi;
00158 }
00159
00160 if (aux_bi_idx<0 && aux_bj_idx<0)
00161 {
00162 if (si<sj)
00163 Columns_data[nbi*block_storage+si+sj*block_size] = val;
00164 else
00165 Columns_data[nbi*block_storage+sj+si*block_size] = val;
00166 }
00167 else if (aux_bi_idx>=0)
00168 column_nbj_bi_data[block_size*sj + si] = val;
00169 else
00170 column_nbi_bj_data[block_size*si + sj] = val;
00171
00172
00173
00174 nonzeros++;
00175 }
00176 }
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 }
00203
00204 SparseGridMtxLL::~SparseGridMtxLL()
00205 {
00206 if (Columns_data) {delete Columns_data; Columns_data = NULL;}
00207 }
00208
00209 void SparseGridMtxLL::AlocateMemoryByPattern(IConectMatrix* bskl)
00210 {
00211 clock_t ts = clock();
00212 this->n_blocks = bskl->N();
00213 columns_data_length = n_blocks*block_storage;
00214 SparseGridColumn* column;
00215
00216 for (long bj=0; bj<n_blocks; bj++)
00217 {
00218 IntArrayList* indexes = bskl->DetachIndexesAboveDiagonalInColumn(bj);
00219 Columns[bj] = column = new SparseGridColumn(indexes);
00220 column->column_start_idx = columns_data_length;
00221 columns_data_length += column->Entries*block_storage;
00222 }
00223
00224 if (columns_data_length==0)
00225 columns_data_length = 1;
00226
00227 long memn = ((long)columns_data_length)*sizeof(double)/1024;
00228 char str[256];
00229
00230 sprintf(str," sparse matrix size : %ld kB..",memn);
00231 Write(str);
00232 this->Columns_data = new double[columns_data_length];
00233
00234 clock_t end = clock();
00235 double duration = (double)(end - ts) / CLOCKS_PER_SEC;
00236 sprintf(str,"%.3f s",duration);
00237 Writeln(str);
00238
00239 }
00240
00241
00242 void SparseGridMtxLL::MultiplyByVector(const LargeVectorAttach& , LargeVectorAttach& )
00243 {
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 }
00266
00267 void SparseGridMtxLL::Factorize()
00268 {
00269
00270 long bi,bj,min_bi_J = 0;
00271
00272
00273 long* p_blockJ_pattern = new long[n_blocks+1];
00274 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00275
00276 long newdotcnount = (long)ceil((double)n_blocks / 24.0);
00277
00278 BlockArith->zero_pivots = 0;
00279
00280 double* cd = this->Columns_data;
00281
00282 {
00283 double* dd = cd;
00284
00285 long Djj = 0;
00286 eMT->act_block = 0;
00287 for (bj=0; bj<n_blocks;bj++)
00288 {
00289 SparseGridColumn &columnJ = *Columns[bj];
00290 long noJentries = columnJ.Entries;
00291 if (noJentries>0)
00292 {
00293 long* columnJentries = columnJ.IndexesUfa->Items;
00294 double* pAkj = cd+columnJ.column_start_idx;
00295 double* pAij = pAkj;
00296
00297
00298 min_bi_J=bj;
00299 for (long i=noJentries-1; i>=0; i--)
00300 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00301
00302
00303 for (long idx_J=0; idx_J<noJentries; idx_J++)
00304 {
00305 bi = columnJentries[idx_J];
00306
00307 SparseGridColumn &columnI = *Columns[bi];
00308 long noIentries = columnI.Entries;
00309
00310 if (noIentries>0)
00311 {
00312 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00313 long* columnIentries = columnI.IndexesUfa->Items;
00314 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00315 {
00316 long idx_K = p_blockJ_pattern[*columnIentry--];
00317 if (idx_K==0) continue;
00318 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00319 }
00320 }
00321 BlockArith->L_BlockSolve(dd+bi*block_storage,pAij);
00322 pAij += block_storage;
00323 }
00324
00325
00326
00327 for (long idx=noJentries-1; idx>=0; idx--)
00328 {
00329 bi = columnJentries[idx];
00330
00331 p_blockJ_pattern[bi] = 0;
00332
00333 int ij = columnJ.column_start_idx + idx*block_storage;
00334 BlockArith->SubATBproduct(dd+Djj,cd+ij,cd+ij);
00335 }
00336 }
00337
00338
00339 BlockArith->LL_Decomposition(dd+Djj);
00340
00341 if ((bj % newdotcnount) == 0)
00342 Write(".");
00343
00344 Djj += block_storage;
00345 eMT->act_block += block_size;
00346 if(eMT->break_flag)
00347 break;
00348 }
00349 }
00350 delete [] p_blockJ_pattern;
00351
00352 ComputeBlocks();
00353 }
00354
00355 void SparseGridMtxLL::Factorize_Incomplete()
00356 {
00357
00358 long bi,bj,min_bi_J = 0;
00359
00360
00361 long* p_blockJ_pattern = new long[n_blocks+1];
00362 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00363
00364 long newdotcnount = (long)ceil((double)n_blocks / 24.0);
00365
00366 BlockArith->zero_pivots = 0;
00367
00368 double* cd = this->Columns_data;
00369
00370 {
00371 double* dd = cd;
00372
00373 long Djj = 0;
00374 eMT->act_block = 0;
00375 for (bj=0; bj<n_blocks;bj++)
00376 {
00377 SparseGridColumn &columnJ = *Columns[bj];
00378 long noJentries = columnJ.Entries;
00379 if (noJentries>0)
00380 {
00381 long* columnJentries = columnJ.IndexesUfa->Items;
00382 double* pAkj = cd+columnJ.column_start_idx;
00383 double* pAij = pAkj;
00384
00385
00386 min_bi_J=bj;
00387 for (long i=noJentries-1; i>=0; i--)
00388 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00389
00390
00391 for (long idx_J=0; idx_J<noJentries; idx_J++)
00392 {
00393 bi = columnJentries[idx_J];
00394
00395 SparseGridColumn &columnI = *Columns[bi];
00396 long noIentries = columnI.Entries;
00397
00398 if (noIentries>0)
00399 {
00400 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00401 long* columnIentries = columnI.IndexesUfa->Items;
00402 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00403 {
00404 long idx_K = p_blockJ_pattern[*columnIentry--];
00405 if (idx_K==0) continue;
00406 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00407 }
00408 }
00409 BlockArith->L_BlockSolve(dd+bi*block_storage,pAij);
00410 pAij += block_storage;
00411 }
00412
00413
00414
00415 for (long idx=noJentries-1; idx>=0; idx--)
00416 {
00417 bi = columnJentries[idx];
00418
00419 p_blockJ_pattern[bi] = 0;
00420
00421 int ij = columnJ.column_start_idx + idx*block_storage;
00422 BlockArith->SubATBproduct(dd+Djj,cd+ij,cd+ij);
00423 }
00424
00425
00426
00427
00428 long maxJentries = columnJ.Entries-1;
00429 if (maxJentries<2)
00430 maxJentries = 2;
00431 if (columnJ.Entries>maxJentries)
00432 {
00433 Array::Copy(columnJ.IndexesUfa->Items+(columnJ.Entries-maxJentries),
00434 columnJ.IndexesUfa->Items,
00435 maxJentries);
00436 columnJ.Entries = maxJentries;
00437 }
00438
00439 }
00440
00441
00442 BlockArith->LL_Decomposition(dd+Djj);
00443
00444 if ((bj % newdotcnount) == 0)
00445 Write(".");
00446
00447 Djj += block_storage;
00448 eMT->act_block += block_size;
00449 if(eMT->break_flag)
00450 break;
00451 }
00452 }
00453 delete [] p_blockJ_pattern;
00454
00455 ComputeBlocks();
00456 }
00457
00458
00459
00460 void SparseGridMtxLL::SchurComplementFactorization(int fixed_blocks)
00461 {
00462
00463 long bi,bj,min_bi_J = 0;
00464
00465 long blocks_to_factor = n_blocks - fixed_blocks;
00466
00467
00468 long* p_blockJ_pattern = new long[n_blocks+1];
00469 memset(p_blockJ_pattern,0,(n_blocks+1)*sizeof(long));
00470
00471
00472
00473 BlockArith->zero_pivots = 0;
00474
00475 double* cd = this->Columns_data;
00476
00477 {
00478 double* dd = cd;
00479
00480 long Djj = 0;
00481 eMT->act_block = 0;
00482 for (bj=0; bj<n_blocks;bj++)
00483 {
00484 SparseGridColumn &columnJ = *Columns[bj];
00485 long noJentries = columnJ.Entries;
00486 if (noJentries>0)
00487 {
00488 long* columnJentries = columnJ.IndexesUfa->Items;
00489 double* pAkj = cd+columnJ.column_start_idx;
00490 double* pAij = pAkj;
00491
00492
00493 min_bi_J=bj;
00494
00495 if (bj<blocks_to_factor)
00496 {
00497 for (long i=noJentries-1; i>=0; i--)
00498 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00499 }
00500 else
00501 {
00502 for (long i=noJentries-1; i>=0; i--)
00503 if (columnJentries[i]<blocks_to_factor)
00504 p_blockJ_pattern[min_bi_J=columnJentries[i]] = ~(i*block_storage);
00505 }
00506
00507
00508 for (long idx_J=0; idx_J<noJentries; idx_J++)
00509 {
00510 bi = columnJentries[idx_J];
00511
00512 SparseGridColumn &columnI = *Columns[bi];
00513 long noIentries = columnI.Entries;
00514
00515 if (noIentries>0)
00516 {
00517 double* pAki = cd+columnI.column_start_idx + (noIentries-1)*block_storage;
00518 long* columnIentries = columnI.IndexesUfa->Items;
00519 for (long* columnIentry = columnIentries+noIentries-1; columnIentry>=columnIentries; pAki-=block_storage)
00520 {
00521 long idx_K = p_blockJ_pattern[*columnIentry--];
00522 if (idx_K==0) continue;
00523 BlockArith->SubATBproduct(pAij,pAki,pAkj+ ~idx_K );
00524 }
00525 }
00526 BlockArith->L_BlockSolve(dd+bi*block_storage,pAij);
00527 pAij += block_storage;
00528 }
00529
00530
00531
00532 for (long idx=noJentries-1; idx>=0; idx--)
00533 {
00534 bi = columnJentries[idx];
00535
00536 p_blockJ_pattern[bi] = 0;
00537
00538 if (bi>=blocks_to_factor)
00539 break;
00540
00541 int ij = columnJ.column_start_idx + idx*block_storage;
00542 BlockArith->SubATBproduct(dd+Djj,cd+ij,cd+ij);
00543 }
00544 }
00545
00546
00547 if (bj<blocks_to_factor)
00548 BlockArith->LL_Decomposition(dd+Djj);
00549
00550
00551
00552
00553 Djj += block_storage;
00554 eMT->act_block += block_size;
00555 if(eMT->break_flag)
00556 break;
00557 }
00558 }
00559 delete [] p_blockJ_pattern;
00560
00561 ComputeBlocks();
00562 }
00563
00564 void SparseGridMtxLL::Solve(double* b, double* x)
00565 {
00566 if (x!=b)
00567 Array::Copy(b,x,n);
00568 SolveLL(x);
00569 }
00570
00571 void SparseGridMtxLL::SolveLV(const LargeVector& b, LargeVector& x)
00572 {
00573 if (&x!=&b)
00574 x.Initialize(b,0,n);
00575 SolveLL(x.DataPtr());
00576 }
00577
00578
00579 void SparseGridMtxLL::ForwardSubstL(double* x,long fixed_blocks)
00580 {
00581 if (this->N()==0) return;
00582 long blocks_to_factor = n_blocks - fixed_blocks;
00583
00584 int bi;
00585 long* ord = this->block_order->order->Items;
00586
00587
00588 for (bi=0; bi<blocks_to_factor; bi++)
00589 {
00590 SparseGridColumn& rowI = *Columns[bi];
00591 int no = rowI.Entries;
00592 if (no>0)
00593 {
00594 double* dst = x+block_size*ord[bi];
00595 double* Aij = Columns_data + rowI.column_start_idx;
00596
00597 long* idxs = rowI.IndexesUfa->Items;
00598
00599 for (int idx = 0; idx <no; idx++,Aij += block_storage)
00600 BlockArith->SubMultTBlockByVector(Aij,x+block_size*ord[*(idxs++)],dst);
00601 }
00602
00603 BlockArith->SubstSolveL(Columns_data+block_storage*bi,x+block_size*ord[bi]);
00604 }
00605 }
00606
00607 void SparseGridMtxLL::BackSubstLT(double* x,long fixed_blocks)
00608 {
00609 if (this->N()==0) return;
00610 long blocks_to_factor = n_blocks - fixed_blocks;
00611
00612 int bi;
00613 long* ord = this->block_order->order->Items;
00614
00615
00616 for (bi=blocks_to_factor-1; bi>=0; bi--)
00617 {
00618 BlockArith->SubstSolveLT(Columns_data+block_storage*bi,x+block_size*ord[bi]);
00619
00620 SparseGridColumn& columnI = *Columns[bi];
00621 int no = columnI.Entries;
00622 if (no>0)
00623 {
00624 double* src = x+block_size*ord[bi];
00625 double* Aij = Columns_data + columnI.column_start_idx;
00626
00627 long* idxs = columnI.IndexesUfa->Items;
00628
00629 for (int idx = 0; idx <no; idx++,Aij += block_storage)
00630 BlockArith->SubMultBlockByVector(Aij,src,x+block_size*ord[*(idxs++)]);
00631 }
00632 }
00633 }
00634
00635
00636 void SparseGridMtxLL::SolveLL(double* x,long fixed_blocks)
00637 {
00638 ForwardSubstL(x,fixed_blocks);
00639 BackSubstLT(x,fixed_blocks);
00640 }
00641
00642 void SparseGridMtxLL::SolveA11(double* x,long fixed_blocks)
00643 {
00644 if (x==0 || fixed_blocks==0)
00645 x = 0;
00646 }
00647
00648 void SparseGridMtxLL::Sub_A21_A11inv(double* x,long fixed_blocks)
00649 {
00650 if (x==0 || fixed_blocks==0)
00651 x = 0;
00652 }
00653
00654 void SparseGridMtxLL::Sub_A11inv_A12(double* x,long fixed_blocks)
00655 {
00656 if (x==0 || fixed_blocks==0)
00657 x = 0;
00658 }
00659
00660 void SparseGridMtxLL::WriteCondensedMatrixA22(double* a,Ordering* mcn,IntArrayList* lncn)
00661 {
00662 long blockSize = BlockSize();
00663 long i,j;
00664 long nbn = lncn->Count;
00665 long aux_bi_idx=0,aux_bj_idx=0;
00666 for (j=0; j<lncn->Count; j++)
00667 {
00668 long nj = mcn->perm->Items[lncn->Items[j]];
00669 long nbj = nj / blockSize;
00670 long nsj = nj % blockSize;
00671 long obj = block_order->perm->Items[nbj];
00672
00673 for (i=0; i<=j; i++)
00674 {
00675 long ni = mcn->perm->Items[lncn->Items[i]];
00676 long nbi = ni / blockSize;
00677 long nsi = ni % blockSize;
00678 long obi = block_order->perm->Items[nbi];
00679
00680 double val = GetValue(obi,obj,nsi,nsj,aux_bi_idx,aux_bj_idx);
00681 a[(j) + (i)*nbn] = val;
00682 a[(i) + (j)*nbn] = val;
00683 }
00684 }
00685 }
00686
00687
00688 DSS_NAMESPASE_END