232 std::vector<int_type> ExtRows;
234 int ExtSize=0,ExtCols=0,MainCols=0;
235 int Nrows=
Matrix_->NumMyRows();
236 int *l2b,*block_offset;
239 int index,col,row_in_block,col_in_block,length,*colind;
242 bool verbose=(bool)(
List_.get(
"output",0) > 0);
255 Matrix_->Comm().ScanSum(&NumBlocks_int_type,&MyMinGID, 1);
257 int_type *MyBlockGIDs=
new int_type[
NumBlocks_];
259 MyBlockGIDs[i]=MyMinGID+i;
271 ExtRows.reserve(ExtSize);
274 if(RowMap.
LID(Blockids_const_ptr<int_type>()[i])==-1)
275 ExtRows.push_back(Blockids_const_ptr<int_type>()[i]);
280 ExtSize=ExtRows.size();
281 RowMap.
Comm().SumAll(&ExtSize,&size_sum,1);
283 if(verbose && !
Matrix_->Comm().MyPID()) printf(
"EpetraExt_PointToBlockDiagPermute: Switching to purely local mode\n");
286 Blockids_ref<int>()[i]=RowMap.
LID(Blockids_const_ptr<int_type>()[i]);
295 block_offset=
new int[Nrows];
300 block_offset[Blockids_const_ptr<int>()[j]]=j-
Blockstart_[i];
301 l2b[Blockids_const_ptr<int>()[j]]=i;
306 for(i=0;i<Nrows;i++) {
307 int block_num = l2b[i];
309 row_in_block=block_offset[i];
310 Matrix_->ExtractMyRowView(i,length,values,colind);
312 for (j = 0; j < length; j++) {
314 if(col < Nrows && l2b[col]==block_num){
316 col_in_block = block_offset[col];
317 index = col_in_block * bsize[block_num] + row_in_block;
318 (*BDMat_)[block_num][index]=values[j];
324 index = row_in_block * bsize[block_num] + row_in_block;
325 (*BDMat_)[block_num][index]=1.0;
331 l2blockid=
new int_type[Nrows];
333 l2blockid[
Blockstart_[l2b[i]]+block_offset[i]]= (int_type)
Matrix_->RowMap().GID64(i);
343 int_type* ExtRowsPtr = ExtRows.size()>0 ? &ExtRows[0] : NULL;
349 TmpMatrix.FillComplete();
351 ExtCols=TmpMatrix.NumMyCols();
356 LocalColIDS=
new int[MainCols+ExtCols];
357 for(i=0;i<MainCols;i++){
358 int_type GID= (int_type)
Matrix_->GCID64(i);
359 int MainLID=RowMap.
LID(GID);
360 if(MainLID!=-1) LocalColIDS[i]=MainLID;
362 int ExtLID=TmpMatrix.LRID(GID);
363 if(ExtLID!=-1) LocalColIDS[i]=Nrows+ExtLID;
364 else LocalColIDS[i]=-1;
369 for(i=0;i<ExtCols;i++){
370 int_type GID= (int_type) TmpMatrix.GCID64(i);
371 int MainLID=RowMap.
LID(GID);
372 if(MainLID!=-1) LocalColIDS[MainCols+i]=MainLID;
374 int ExtLID=TmpMatrix.LRID(GID);
375 if(ExtLID!=-1) LocalColIDS[MainCols+i]=Nrows+ExtLID;
376 else LocalColIDS[MainCols+i]=-1;
381 l2b=
new int[Nrows+ExtSize];
382 block_offset=
new int[Nrows+ExtSize];
387 for(i=0;i<Nrows+ExtSize;i++) block_offset[i]=l2b[i]=-1;
391 int LID=RowMap.
LID(Blockids_const_ptr<int_type>()[j]);
392 if(LID==-1) {LID=Nrows+ext_idx;ext_idx++;}
399 for(i=0;i<Nrows;i++) {
400 int block_num = l2b[i];
402 row_in_block=block_offset[i];
403 Matrix_->ExtractMyRowView(i,length,values,colind);
405 for (j = 0; j < length; j++) {
406 col = LocalColIDS[colind[j]];
407 if(col!=-1 && l2b[col]==block_num){
409 col_in_block = block_offset[col];
410 index = col_in_block * bsize[block_num] + row_in_block;
411 (*BDMat_)[block_num][index]=values[j];
417 index = row_in_block * bsize[block_num] + row_in_block;
418 (*BDMat_)[block_num][index]=values[j];
424 for(i=0;i<ExtSize;i++) {
425 int block_num = l2b[Nrows+i];
427 row_in_block=block_offset[Nrows+i];
428 TmpMatrix.ExtractMyRowView(i,length,values,colind);
430 for (j = 0; j < length; j++) {
431 col = LocalColIDS[MainCols+colind[j]];
432 if(col!=-1 && l2b[col]==block_num){
434 col_in_block = block_offset[col];
435 index = col_in_block * bsize[block_num] + row_in_block;
436 (*BDMat_)[block_num][index]=values[j];
442 index = row_in_block * bsize[block_num] + row_in_block;
443 (*BDMat_)[block_num][index]=1.0;
450 for(i=0;i<Nrows;i++){
452 if(bkid>-1) l2blockid[
Blockstart_[bkid]+block_offset[i]]= (int_type) RowMap.
GID64(i);
455 for(i=0;i<ExtSize;i++)
456 l2blockid[
Blockstart_[l2b[Nrows+i]]+block_offset[Nrows+i]]= (int_type) TmpMatrix.GRID64(i);
464 Teuchos::ParameterList dummy,inList;
465 inList=
List_.get(
"blockdiagmatrix: list",dummy);
476 delete [] LocalColIDS;
477 delete [] block_offset;
481 delete [] MyBlockGIDs;