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