Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_VbrMatrix.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42#include "Epetra_ConfigDefs.h"
43#include "Epetra_VbrMatrix.h"
44#include "Epetra_BlockMap.h"
45#include "Epetra_Map.h"
46#include "Epetra_Import.h"
47#include "Epetra_Export.h"
48#include "Epetra_Vector.h"
49#include "Epetra_MultiVector.h"
50#include "Epetra_Comm.h"
51#include "Epetra_Distributor.h"
53
54#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
55// FIXME long long : whole file
56
57//==============================================================================
58Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& rowMap, int *NumBlockEntriesPerRow)
59 : Epetra_DistObject(rowMap, "Epetra::VbrMatrix"),
62 Graph_(0),
63 Allocated_(false),
64 StaticGraph_(false),
65 constructedWithFilledGraph_(false),
66 matrixFillCompleteCalled_(false),
67 NumMyBlockRows_(rowMap.NumMyElements()),
68 CV_(CV),
69 NormInf_(0.0),
70 NormOne_(0.0),
71 NormFrob_(0.0),
72 squareFillCompleteCalled_(false)
73{
75 Graph_ = new Epetra_CrsGraph(CV, rowMap, NumBlockEntriesPerRow);
76
77#ifdef NDEBUG
78 (void) Allocate();
79#else
80 // Don't declare 'err' unless we actually use it (in the assert(),
81 // which gets defined away in a release build).
82 int err = Allocate();
83 assert( err == 0 );
84#endif // NDEBUG
85}
86
87//==============================================================================
88Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& rowMap, int NumBlockEntriesPerRow)
89 : Epetra_DistObject(rowMap, "Epetra::VbrMatrix"),
92 Graph_(0),
93 Allocated_(false),
94 StaticGraph_(false),
95 constructedWithFilledGraph_(false),
96 matrixFillCompleteCalled_(false),
97 NumMyBlockRows_(rowMap.NumMyElements()),
98 CV_(CV),
99 squareFillCompleteCalled_(false)
100{
102 Graph_ = new Epetra_CrsGraph(CV, rowMap, NumBlockEntriesPerRow);
103
104#ifdef NDEBUG
105 (void) Allocate();
106#else
107 // Don't declare 'err' unless we actually use it (in the assert(),
108 // which gets defined away in a release build).
109 int err = Allocate();
110 assert( err == 0 );
111#endif // NDEBUG
112}
113//==============================================================================
115 const Epetra_BlockMap& colMap, int *NumBlockEntriesPerRow)
116 : Epetra_DistObject(rowMap, "Epetra::VbrMatrix"),
118 Epetra_BLAS(),
119 Graph_(0),
120 Allocated_(false),
121 StaticGraph_(false),
122 constructedWithFilledGraph_(false),
123 matrixFillCompleteCalled_(false),
124 NumMyBlockRows_(rowMap.NumMyElements()),
125 CV_(CV),
126 squareFillCompleteCalled_(false)
127{
129 Graph_ = new Epetra_CrsGraph(CV, rowMap, colMap, NumBlockEntriesPerRow);
130
131#ifdef NDEBUG
132 (void) Allocate();
133#else
134 // Don't declare 'err' unless we actually use it (in the assert(),
135 // which gets defined away in a release build).
136 int err = Allocate();
137 assert( err == 0 );
138#endif // NDEBUG
139}
140
141//==============================================================================
143 const Epetra_BlockMap& colMap, int NumBlockEntriesPerRow)
144 : Epetra_DistObject(rowMap, "Epetra::VbrMatrix"),
146 Epetra_BLAS(),
147 Graph_(0),
148 Allocated_(false),
149 StaticGraph_(false),
150 constructedWithFilledGraph_(false),
151 matrixFillCompleteCalled_(false),
152 NumMyBlockRows_(rowMap.NumMyElements()),
153 CV_(CV),
154 squareFillCompleteCalled_(false)
155{
157 Graph_ = new Epetra_CrsGraph(CV, rowMap, colMap, NumBlockEntriesPerRow);
158
159#ifdef NDEBUG
160 (void) Allocate();
161#else
162 // Don't declare 'err' unless we actually use it (in the assert(),
163 // which gets defined away in a release build).
164 int err = Allocate();
165 assert( err == 0 );
166#endif // NDEBUG
167}
168
169//==============================================================================
171 : Epetra_DistObject(graph.RowMap(), "Epetra::VbrMatrix"),
173 Epetra_BLAS(),
174 Graph_(new Epetra_CrsGraph(graph)),
175 Allocated_(false),
176 StaticGraph_(true),
177 constructedWithFilledGraph_(false),
178 matrixFillCompleteCalled_(false),
179 NumMyBlockRows_(graph.RowMap().NumMyElements()),
180 CV_(CV),
181 squareFillCompleteCalled_(false)
182{
185
186#ifdef NDEBUG
187 (void) Allocate();
188#else
189 // Don't declare 'err' unless we actually use it (in the assert(),
190 // which gets defined away in a release build).
191 int err = Allocate();
192 assert( err == 0 );
193#endif // NDEBUG
194}
195
196//==============================================================================
198 : Epetra_DistObject(Source),
200 Epetra_BLAS(),
201 Graph_(new Epetra_CrsGraph(Source.Graph())),
202 Allocated_(Source.Allocated_),
203 StaticGraph_(true),
204 UseTranspose_(Source.UseTranspose_),
205 constructedWithFilledGraph_(Source.constructedWithFilledGraph_),
206 matrixFillCompleteCalled_(Source.matrixFillCompleteCalled_),
207 NumMyBlockRows_(0),
208 CV_(Copy),
209 HavePointObjects_(false),
210 squareFillCompleteCalled_(false)
211 {
212
214 operator=(Source);
215}
216
217//==============================================================================
219{
220 if (this == &src) {
221 return(*this);
222 }
223
224 DeleteMemory();
225
230 CV_ = src.CV_;
231
233
234 //our Graph_ member was delete in DeleteMemory() so we don't need to
235 //delete it now before re-creating it.
236 Graph_ = new Epetra_CrsGraph(src.Graph());
237
238#ifdef NDEBUG
239 (void) Allocate();
240#else
241 // Don't declare 'err' unless we actually use it (in the assert(),
242 // which gets defined away in a release build).
243 int err = Allocate();
244 assert( err == 0 );
245#endif // NDEBUG
246
247 int i, j;
248
249#if 0
250 bool ConstantShape = true;
251 const int NOTSETYET = -13 ;
252 int MyLda = NOTSETYET ;
253 int MyColDim = NOTSETYET ;
254 int MyRowDim = NOTSETYET ;
255 // int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
256 // if (ierr) EPETRA_CHK_ERR(ierr);
257 // if ( ierr != 0 ) ConstantShape = false ; // Do we really need this?
258 assert( ConstantShape ) ;
259 for (i=0; i<NumMyBlockRows_; i++) {
260 int NumBlockEntries = NumBlockEntriesPerRow_[i];
261 for (j=0; j < NumBlockEntries; j++) {
262 Epetra_SerialDenseMatrix* ThisBlock = src.Entries_[i][j];
263 if ( MyLda == NOTSETYET ) {
264 MyLda = ThisBlock->LDA() ;
265 MyColDim = ThisBlock->ColDim() ;
266 MyRowDim = ThisBlock->RowDim() ;
267 } else {
268 if ( MyLda != ThisBlock->LDA() ) ConstantShape = false ;
269 if ( MyRowDim != ThisBlock->RowDim() ) ConstantShape = false ;
270 if ( MyColDim != ThisBlock->ColDim() ) ConstantShape = false ;
271 }
272 }
273 }
274 if ( false && ConstantShape ) {
275
277
278 All_Values_ = new double[NumMyNonzeros];
280 for (i=0; i<NumMyBlockRows_; i++) {
281 double* Values_ThisBlockRow = All_Values_ ;
282 int NumBlockEntries = NumBlockEntriesPerRow_[i];
283 for (j=0; j < NumBlockEntries; j++) {
284 double* Values_ThisBlockEntry = All_Values_ ;
285 Epetra_SerialDenseMatrix* M_SDM = src.Entries_[i][j] ;
286 for ( int kk = 0 ; kk < MyColDim ; kk++ ) {
287 for ( int ll = 0 ; ll < MyRowDim ; ll++ ) {
288 *All_Values_ = (*M_SDM)(ll,kk);
289 All_Values_++ ;
290 }
291 }
292 Entries_[i][j] = new Epetra_SerialDenseMatrix(View, Values_ThisBlockEntry,
293 MyLda, MyRowDim, MyColDim );
294 }
295 }
296 } else {
297#endif
298 for (i=0; i<NumMyBlockRows_; i++) {
299 int NumBlockEntries = NumBlockEntriesPerRow_[i];
300 for (j=0; j < NumBlockEntries; j++) {
301 //if src.Entries_[i][j] != 0, set Entries_[i][j] to be
302 //a copy of it.
303 Entries_[i][j] = src.Entries_[i][j] != 0 ?
304 new Epetra_SerialDenseMatrix(*(src.Entries_[i][j])) : 0;
305 }
306#if 0
307 }
308#endif
309 }
310
311 if ( src.StorageOptimized() ) this->OptimizeStorage() ;
312
313 return( *this );
314}
315
316//==============================================================================
317void Epetra_VbrMatrix::InitializeDefaults() { // Initialize all attributes that have trivial default values
318
319 UseTranspose_ = false;
320 Entries_ = 0;
321 All_Values_ = 0;
323 NormInf_ = -1.0;
324 NormOne_ = -1.0;
325 NormFrob_ = -1.0;
326 ImportVector_ = 0;
327
330 Indices_ = 0;
333
334 // State variables needed for constructing matrix entry-by-entry
335
336 TempRowDims_ = 0;
337 TempEntries_ = 0;
338 LenTemps_ = 0;
339 CurBlockRow_ = 0;
342 CurEntry_ = -1; // Set to -1 to allow a simple sanity check when submitting entries
343 CurIndicesAreLocal_ = false;
345
346 // State variables needed for extracting entries
348 CurExtractEntry_ = -1; // Set to -1 to allow a simple sanity check when extracting entries
351 CurExtractView_ = false;
352 CurRowDim_ = 0;
353
354 // State variable for extracting block diagonal entries
355 CurBlockDiag_ = -1; // Set to -1 to allow a simple sanity check when extracting entries
356
357 // Attributes that support the Epetra_RowMatrix and Epetra_Operator interfaces
363 HavePointObjects_ = false;
364 StorageOptimized_ = false ;
365
366 OperatorX_ = 0;
367 OperatorY_ = 0;
368
369 return;
370}
371
372//==============================================================================
374
375 int i, j;
376
377 // Set direct access pointers to graph info (needed for speed)
381
383
385
386
387 // Allocate Entries array
389 // Allocate and initialize entries
390 for (i=0; i<NumMyBlockRows_; i++) {
391 int NumAllocatedBlockEntries = NumAllocatedBlockEntriesPerRow_[i];
392
393 if (NumAllocatedBlockEntries > 0) {
394 Entries_[i] = new Epetra_SerialDenseMatrix*[NumAllocatedBlockEntries];
395 for (j=0; j < NumAllocatedBlockEntries; j++) {
396 Entries_[i][j] = 0;
397 }
398 }
399 else {
400 Entries_[i] = 0;
401 }
402 }
403 SetAllocated(true);
404 return(0);
405}
406//==============================================================================
408{
409 DeleteMemory();
410}
411
412//==============================================================================
414{
415 int i;
417 for (i=0; i<NumMyBlockRows_; i++) {
418 int NumAllocatedBlockEntries = NumAllocatedBlockEntriesPerRow_[i];
419
420 if (NumAllocatedBlockEntries >0) {
421
422 for (int j=0; j < NumAllocatedBlockEntries; j++) {
423 if (Entries_[i][j]!=0) {
424 delete Entries_[i][j];
425 }
426 }
427 delete [] Entries_[i];
428 }
429 }
430
431 if (All_Values_Orig_!=0) delete [] All_Values_Orig_;
433
434 if (Entries_!=0) delete [] Entries_;
435 Entries_ = 0;
436
437 if (ImportVector_!=0) delete ImportVector_;
438 ImportVector_ = 0;
439
440 NumMyBlockRows_ = 0;
441
442 if (LenTemps_>0) {
443 delete [] TempRowDims_;
444 delete [] TempEntries_;
445 }
446
447 // Delete any objects related to supporting the RowMatrix and Operator interfaces
448 if (HavePointObjects_) {
449#if 1
453#else
454 // this can fail, see bug # 1114
455 if (!RowMatrixColMap().SameAs(RowMatrixRowMap())) delete RowMatrixColMap_;
457 if (!OperatorRangeMap().SameAs(RowMatrixRowMap())) delete OperatorRangeMap_;
458#endif
459 delete RowMatrixRowMap_;
460 delete RowMatrixImporter_;
461 HavePointObjects_ = false;
462 }
463
464 if (OperatorX_!=0) {
465 delete OperatorX_;
466 delete OperatorY_;
467 }
468
469 InitializeDefaults(); // Reset all basic pointers to zero
470 Allocated_ = false;
471
472 delete Graph_; // We created the graph, so must delete it.
473 Graph_ = 0;
474}
475
476//==============================================================================
477int Epetra_VbrMatrix::PutScalar(double ScalarConstant)
478{
479 if (!Allocated_) return(0);
480
481 for (int i=0; i<NumMyBlockRows_; i++) {
482 int NumBlockEntries = NumBlockEntriesPerRow_[i];
483 int RowDim = ElementSizeList_[i];
484 for (int j=0; j< NumBlockEntries; j++) {
485 int LDA = Entries_[i][j]->LDA();
486 int ColDim = Entries_[i][j]->N();
487 for (int col=0; col < ColDim; col++) {
488 double * Entries = Entries_[i][j]->A()+col*LDA;
489 for (int row=0; row < RowDim; row++)
490 *Entries++ = ScalarConstant;
491 }
492 }
493 }
494 NormOne_ = -1.0; // Reset Norm so it will be recomputed.
495 NormInf_ = -1.0; // Reset Norm so it will be recomputed.
496 NormFrob_ = -1.0;
497 return(0);
498}
499
500//==============================================================================
501int Epetra_VbrMatrix::Scale(double ScalarConstant)
502{
503 for (int i=0; i<NumMyBlockRows_; i++) {
504 int NumBlockEntries = NumBlockEntriesPerRow_[i];
505 int RowDim = ElementSizeList_[i];
506 for (int j=0; j< NumBlockEntries; j++) {
507 int LDA = Entries_[i][j]->LDA();
508 int ColDim = Entries_[i][j]->N();
509 for (int col=0; col < ColDim; col++) {
510 double * Entries = Entries_[i][j]->A()+col*LDA;
511 for (int row=0; row < RowDim; row++)
512 *Entries++ *= ScalarConstant;
513 }
514 }
515 }
516 NormOne_ = -1.0; // Reset Norm so it will be recomputed.
517 NormInf_ = -1.0; // Reset Norm so it will be recomputed.
518 NormFrob_ = -1.0;
519 return(0);
520}
521//==========================================================================
522int Epetra_VbrMatrix::BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int * BlockIndices) {
523
524 if (IndicesAreLocal()) EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
526 int LocalBlockRow = LRID(BlockRow); // Find local row number for this global row index
527
528 bool indicesAreLocal = false;
529 EPETRA_CHK_ERR(BeginInsertValues(LocalBlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
530 return(0);
531}
532
533//==========================================================================
534int Epetra_VbrMatrix::BeginInsertMyValues(int BlockRow, int NumBlockEntries,
535 int * BlockIndices)
536{
537 if (IndicesAreGlobal()) EPETRA_CHK_ERR(-2); // Cannot insert local values until MakeIndicesLocal() is called either directly or through FillComplete()
539 bool indicesAreLocal = true;
540 EPETRA_CHK_ERR(BeginInsertValues(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
541 return(0);
542
543}
544
545//==========================================================================
546int Epetra_VbrMatrix::BeginInsertValues(int BlockRow, int NumBlockEntries,
547 int * BlockIndices, bool indicesAreLocal)
548{
549 if (StaticGraph()) EPETRA_CHK_ERR(-2); // If the matrix graph is fully constructed, we cannot insert new values
550
551 int ierr = 0;
552
553 if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) {
554 EPETRA_CHK_ERR(-1); // Not in BlockRow range
555 }
556 if (CV_==View && Entries_[BlockRow]!=0) ierr = 2; // This row has be defined already. Issue warning.
557 if (IndicesAreContiguous()) EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
558
559 // Set up pointers, make sure enough temp space for this round of submits
560
561 Epetra_CombineMode SubmitMode = Insert;
562
563 EPETRA_CHK_ERR(ierr);
564 EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal, SubmitMode));
565 return(0);
566}
567//==========================================================================
568int Epetra_VbrMatrix::BeginReplaceGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
569
570 BlockRow = LRID(BlockRow); // Normalize row range
571 bool indicesAreLocal = false;
572 EPETRA_CHK_ERR(BeginReplaceValues(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
573 return(0);
574}
575
576//==========================================================================
577int Epetra_VbrMatrix::BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
578
580
581 bool indicesAreLocal = true;
582 EPETRA_CHK_ERR(BeginReplaceValues(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
583 return(0);
584}
585
586//==========================================================================
588 int NumBlockEntries,
589 int *BlockIndices,
590 bool indicesAreLocal)
591{
592 if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) EPETRA_CHK_ERR(-1); // Not in BlockRow range
593
594 Epetra_CombineMode SubmitMode = Zero; // This is a misuse of Zero mode, fix it later
595 EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal, SubmitMode));
596 return(0);
597}
598
599//==========================================================================
600int Epetra_VbrMatrix::BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
601
602 BlockRow = LRID(BlockRow); // Normalize row range
603 bool indicesAreLocal = false;
604 EPETRA_CHK_ERR(BeginSumIntoValues(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
605 return(0);
606}
607
608//==========================================================================
609int Epetra_VbrMatrix::BeginSumIntoMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
610
611 bool indicesAreLocal = true;
612 EPETRA_CHK_ERR(BeginSumIntoValues(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal));
613 return(0);
614}
615
616//==========================================================================
617int Epetra_VbrMatrix::BeginSumIntoValues(int BlockRow, int NumBlockEntries,
618 int *BlockIndices, bool indicesAreLocal) {
619
620 if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) EPETRA_CHK_ERR(-1); // Not in BlockRow range
621
622 Epetra_CombineMode SubmitMode = Add;
623 EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, indicesAreLocal, SubmitMode));
624 return(0);
625}
626
627//==========================================================================
628int Epetra_VbrMatrix::SetupForSubmits(int BlockRow, int NumBlockEntries,
629 int * BlockIndices,
630 bool indicesAreLocal,
631 Epetra_CombineMode SubmitMode) {
632
633 if (NumBlockEntries>LenTemps_) {
634 if (LenTemps_>0) {
635 delete [] TempRowDims_;
636 delete [] TempEntries_;
637 }
638 TempRowDims_ = new int[NumBlockEntries];
639 TempEntries_ = new Epetra_SerialDenseMatrix*[NumBlockEntries];
640 LenTemps_ = NumBlockEntries;
641 }
642
643 CurBlockRow_ = BlockRow;
644 CurNumBlockEntries_ = NumBlockEntries;
645 CurBlockIndices_ = BlockIndices;
646 CurEntry_ = 0;
647 CurIndicesAreLocal_ = indicesAreLocal;
648 CurSubmitMode_ = SubmitMode;
649
650 return(0);
651}
652
653//==========================================================================
654int Epetra_VbrMatrix::DirectSubmitBlockEntry(int GlobalBlockRow, int GlobalBlockCol,
655 const double *values, int LDA,
656 int NumRows, int NumCols, bool sum_into)
657{
658 int ierr = 0;
659 int LocalBlockRow = LRID(GlobalBlockRow); // Normalize row range
660 int Loc;
661 bool found = Graph_->FindGlobalIndexLoc(LocalBlockRow, GlobalBlockCol,0,Loc);
662 if (found) {
663 if (Entries_[LocalBlockRow][Loc]==0) {
664 Entries_[LocalBlockRow][Loc] =
665 new Epetra_SerialDenseMatrix(Copy, const_cast<double*>(values), LDA, NumRows, NumCols, false);
666 }
667 else {
668 Epetra_SerialDenseMatrix* target = Entries_[LocalBlockRow][Loc];
669
670 target->CopyMat(values, LDA, NumRows, NumCols,
671 target->A_, target->LDA_, sum_into);
672 }
673 }
674 else {
675 ierr = -2;
676 }
677 EPETRA_CHK_ERR(ierr);
678 return ierr;
679}
680
681//==========================================================================
682int Epetra_VbrMatrix::SubmitBlockEntry(double *values, int LDA,
683 int NumRows, int NumCols)
684{
685 if (CurEntry_==-1) EPETRA_CHK_ERR(-1); // This means that a Begin routine was not called
686 if (CurEntry_>=CurNumBlockEntries_) EPETRA_CHK_ERR(-4); // Exceeded the number of entries that can be submitted
687
688 // Fill up temp space with entry
689
690 TempRowDims_[CurEntry_] = NumRows;
692 NumRows, NumCols, false);
693 CurEntry_++;
694
695 return(0);
696}
697
698//==========================================================================
700{
701 return SubmitBlockEntry( Mat.A(), Mat.LDA(), Mat.M(), Mat.N() );
702}
703
704//==========================================================================
706
707 if (CurEntry_!=CurNumBlockEntries_) EPETRA_CHK_ERR(-6); // Did not submit right number of entries
708
709 if (CurSubmitMode_==Insert) {
711 }
712 else {
714 }
715 NormOne_ = -1.0; // Reset Norm so it will be recomputed.
716 NormInf_ = -1.0; // Reset Norm so it will be recomputed.
717 NormFrob_ = -1.0;
718 return(0);
719}
720
721//==========================================================================
723{
724 int j;
725 int ierr = 0;
726 int Loc;
727
728 int RowDim = ElementSizeList_[CurBlockRow_];
729
730 bool SumInto = (CurSubmitMode_==Add);
731
732 for (j=0; j<CurNumBlockEntries_; j++) {
733 int BlockIndex = CurBlockIndices_[j];
734
735 bool code = false;
737 code =Graph_->FindMyIndexLoc(CurBlockRow_,BlockIndex,j,Loc);
738 }
739 else {
740 code =Graph_->FindGlobalIndexLoc(CurBlockRow_,BlockIndex,j,Loc);
741 }
742
743 bool need_to_delete_temp_entry = true;
744
745 if (code) {
747
748 int ColDim = src->N();
749
750 if (Entries_[CurBlockRow_][Loc]==0) {
751 Entries_[CurBlockRow_][Loc] = src;
752 need_to_delete_temp_entry = false;
753 }
754 else {
756
757 target->CopyMat(src->A_, src->LDA_, RowDim, ColDim,
758 target->A_, target->LDA_, SumInto);
759 }
760 }
761 else {
762 ierr=2; // Block Discarded, Not Found
763 }
764
765 if (need_to_delete_temp_entry) {
766 delete TempEntries_[j];
767 }
768 }
769
770 EPETRA_CHK_ERR(ierr);
771 return(0);
772}
773
774//==========================================================================
776{
777 int ierr = 0;
778 int j;
779
780 int NumValidBlockIndices = CurNumBlockEntries_;
781 int * ValidBlockIndices = new int[CurNumBlockEntries_];
782 for ( j=0; j < CurNumBlockEntries_; ++j ) ValidBlockIndices[j] = j;
783
784 if( Graph_->HaveColMap() ) { //test and discard indices not in ColMap
785 NumValidBlockIndices = 0;
786 const Epetra_BlockMap& map = Graph_->ColMap();
787
788 for ( j = 0; j < CurNumBlockEntries_; ++j ) {
789 bool myID = CurIndicesAreLocal_ ?
791
792 if( myID ) {
793 ValidBlockIndices[ NumValidBlockIndices++ ] = j;
794 }
795 else ierr=2; // Discarding a Block not found in ColMap
796 }
797 }
798
799 int oldNumBlocks = NumBlockEntriesPerRow_[CurBlockRow_];
800 int oldNumAllocBlocks = NumAllocatedBlockEntriesPerRow_[CurBlockRow_];
801
802 // Update graph
804
805 int newNumAllocBlocks = NumAllocatedBlockEntriesPerRow_[CurBlockRow_];
806
807 if (newNumAllocBlocks > oldNumAllocBlocks) {
808 ierr = 3;//positive warning code indicates we expanded allocated memory
809 Epetra_SerialDenseMatrix** tmp_Entries = new Epetra_SerialDenseMatrix*[newNumAllocBlocks];
810 for(j=0; j<oldNumBlocks; ++j) tmp_Entries[j] = Entries_[CurBlockRow_][j];
811 for(j=oldNumBlocks; j<newNumAllocBlocks; ++j) tmp_Entries[j] = NULL;
812 delete [] Entries_[CurBlockRow_];
813 Entries_[CurBlockRow_] = tmp_Entries;
814 }
815
816 int start = oldNumBlocks;
817 int stop = start + NumValidBlockIndices;
818 int NumAllocatedEntries = NumAllocatedBlockEntriesPerRow_[CurBlockRow_];
819
820 if (stop <= NumAllocatedEntries) {
821 for (j=start; j<stop; j++) {
823 *(TempEntries_[ValidBlockIndices[j-start]]);
824
826 mat.LDA(),
827 mat.M(),
828 mat.N());
829 }
830 }
831 else {
832 ierr = -4;
833 }
834
835
836 delete [] ValidBlockIndices;
837
838 for (j=0; j<CurNumBlockEntries_; ++j) {
839 delete TempEntries_[j];
840 }
841
842 EPETRA_CHK_ERR(ierr);
843
844 return(0);
845}
846
847//=============================================================================
848int Epetra_VbrMatrix::CopyMat(double * A, int LDA, int NumRows, int NumCols,
849 double * B, int LDB, bool SumInto) const
850{
851 int i, j;
852 double * ptr1 = B;
853 double * ptr2;
854
855 if (LDB<NumRows) EPETRA_CHK_ERR(-1); // Stride of B is not large enough
856 if (LDA<NumRows) EPETRA_CHK_ERR(-1); // Stride of A is not large enough
857
858 if (SumInto) { // Add to existing values
859 for (j=0; j<NumCols; j++) {
860 ptr1 = B + j*LDB;
861 ptr2 = A + j*LDA;
862 for (i=0; i<NumRows; i++) *ptr1++ += *ptr2++;
863 }
864 }
865 else { // Replace values
866 for (j=0; j<NumCols; j++) {
867 ptr1 = B + j*LDB;
868 ptr2 = A + j*LDA;
869 for (i=0; i<NumRows; i++) *ptr1++ = *ptr2++;
870 }
871 }
872 return(0);
873}
874//==========================================================================
878 return(0);
879}
880
881//==========================================================================
883 const Epetra_BlockMap& range_map)
884{
885 int returnValue = 0;
886
887 if (Graph_->Filled()) {
889 returnValue = 2;
890 }
891 }
892
893 if(!StaticGraph()) {
894 EPETRA_CHK_ERR(Graph_->MakeIndicesLocal(domain_map, range_map));
895 }
896
897 EPETRA_CHK_ERR(SortEntries()); // Sort column entries from smallest to largest
898 EPETRA_CHK_ERR(MergeRedundantEntries()); // Get rid of any redundant index values
899
900 if(!StaticGraph()) {
901 EPETRA_CHK_ERR(Graph_->FillComplete(domain_map, range_map));
902 }
903
904 // NumMyCols_ = Graph_->NumMyCols(); // Redefine based on local number of cols
905
907
909 if (DomainMap().NumGlobalElements64() != RangeMap().NumGlobalElements64()) {
910 returnValue = 3;
911 }
913 EPETRA_CHK_ERR(returnValue);
914 }
915
916 return(returnValue);
917}
918
919//==========================================================================
922 return(0);
923}
924
925//==========================================================================
927 EPETRA_CHK_ERR(FillComplete(*domainMap, *rangeMap));
928 return(0);
929}
930
931//==========================================================================
933
935 if (Sorted()) return(0);
936
937 // For each row, sort column entries from smallest to largest.
938 // Use shell sort. Stable sort so it is fast if indices are already sorted.
939
940
941 for (int i=0; i<NumMyBlockRows_; i++){
942
943 Epetra_SerialDenseMatrix ** Entries = Entries_[i];
944 int NumEntries = NumBlockEntriesPerRow_[i];
945 int * Indices = Indices_[i];
946 int n = NumEntries;
947 int m = n/2;
948
949 while (m > 0) {
950 int max = n - m;
951 for (int j=0; j<max; j++)
952 {
953 for (int k=j; k>=0; k-=m)
954 {
955 if (Indices[k+m] >= Indices[k])
956 break;
957 Epetra_SerialDenseMatrix *dtemp = Entries[k+m];
958 Entries[k+m] = Entries[k];
959 Entries[k] = dtemp;
960
961 int itemp = Indices[k+m];
962 Indices[k+m] = Indices[k];
963 Indices[k] = itemp;
964 }
965 }
966 m = m/2;
967 }
968 }
969 Graph_->SetSorted(true); // This also sorted the graph
970 return(0);
971}
972
973//==========================================================================
975{
976
977 if (NoRedundancies()) return(0);
978 if (!Sorted()) EPETRA_CHK_ERR(-1); // Must have sorted entries
979
980 // For each row, remove column indices that are repeated.
981 // Also, determine if matrix is upper or lower triangular or has no diagonal
982 // Note: This function assumes that SortEntries was already called.
983
984 bool SumInto = true;
985 for (int i=0; i<NumMyBlockRows_; i++){
986 int NumEntries = NumBlockEntriesPerRow_[i];
987 if (NumEntries>1) {
988 Epetra_SerialDenseMatrix ** const Entries = Entries_[i];
989 int * const Indices = Indices_[i];
990 int RowDim = ElementSizeList_[i];
991 int curEntry =0;
992 Epetra_SerialDenseMatrix* curBlkEntry = Entries[0];
993 for (int k=1; k<NumEntries; k++) {
994 if (Indices[k]==Indices[k-1]) {
995 if (curBlkEntry->M() != Entries[curEntry]->M() ||
996 curBlkEntry->N() != Entries[curEntry]->N() ||
997 curBlkEntry->LDA() != Entries[curEntry]->LDA()) {
998 std::cerr << "Epetra_VbrMatrix ERROR, two dense-matrix contributions to the same column-index have different sizes: ("<<curBlkEntry->M()<<"x"<<curBlkEntry->N()<<") and ("<<Entries[curEntry]->M()<<"x"<<Entries[curEntry]->N()<<")"<<std::endl;
999 EPETRA_CHK_ERR(-1);
1000 }
1001
1002 CopyMat(Entries[k]->A(), Entries[k]->LDA(), RowDim, Entries[k]->N(),
1003 curBlkEntry->A(), curBlkEntry->LDA(), SumInto);
1004 }
1005 else {
1006 curBlkEntry = Entries[++curEntry];
1007 if (curEntry!=k) {
1008 curBlkEntry->Shape(RowDim,Entries[k]->N());
1009 EPETRA_CHK_ERR(CopyMat(Entries[k]->A(), Entries[k]->LDA(), RowDim, Entries[k]->N(),
1010 curBlkEntry->A(), curBlkEntry->LDA(), false));
1011 }
1012 }
1013 }
1014 }
1015 }
1016
1017 EPETRA_CHK_ERR(Graph_->RemoveRedundantIndices()); // Remove redundant indices and then return
1018 return(0);
1019}
1020
1021//==========================================================================
1023
1024 if (StorageOptimized()) return(0); // Have we been here before?
1025
1026 // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << std::endl ;
1027
1028 bool ConstantShape = true;
1029 int i,j;
1030 const int NOTSETYET = -13 ;
1031 int MyLda = NOTSETYET ;
1032 int MyColDim = NOTSETYET ;
1033 int MyRowDim = NOTSETYET ;
1034 //
1035 // I don't know why the underlying graph musthave optimized storage, but
1036 // the test for optimized storage depends upon the graph hainvg optimized
1037 // storage and that is enough for me.
1038 //
1039 // Ideally, we would like to have optimized storage for the graph. But,
1040 // for now, this is commented out until bug 1097 is fixed
1041 //
1042 // int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
1043 // if (ierr) EPETRA_CHK_ERR(ierr);
1044 // if ( ierr != 0 ) ConstantShape = false ;
1045 if( ConstantShape ) {
1046 for (i=0; i<NumMyBlockRows_; i++) {
1047 int NumBlockEntries = NumBlockEntriesPerRow_[i];
1048 for (j=0; j < NumBlockEntries; j++) {
1049 Epetra_SerialDenseMatrix* ThisBlock = Entries_[i][j];
1050 if ( MyLda == NOTSETYET ) {
1051 MyLda = ThisBlock->LDA() ;
1052 MyColDim = ThisBlock->ColDim() ;
1053 MyRowDim = ThisBlock->RowDim() ;
1054 } else {
1055 if ( MyLda != ThisBlock->LDA() ) ConstantShape = false ;
1056 if ( MyRowDim != ThisBlock->RowDim() ) ConstantShape = false ;
1057 if ( MyColDim != ThisBlock->ColDim() ) ConstantShape = false ;
1058 }
1059 }
1060 }
1061 }
1062 // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << std::endl ;
1063
1064 if ( ConstantShape ) {
1065
1066 int numMyNonzeros = Graph_->NumMyEntries();
1067 int coef_len = MyColDim*MyRowDim*numMyNonzeros;
1068 All_Values_ = new double[coef_len];
1070 for (i=0; i<NumMyBlockRows_; i++) {
1071 int NumBlockEntries = NumBlockEntriesPerRow_[i];
1072 for (j=0; j < NumBlockEntries; j++) {
1073 double* Values_ThisBlockEntry = All_Values_ ;
1074 Epetra_SerialDenseMatrix* M_SDM = Entries_[i][j] ;
1075 for ( int kk = 0 ; kk < MyColDim ; kk++ ) {
1076 for ( int ll = 0 ; ll < MyRowDim ; ll++ ) {
1077 *All_Values_ = (*M_SDM)(ll,kk) ;
1078 All_Values_++ ;
1079 }
1080 }
1081 // Entries_[i][j] = new Epetra_SerialDenseMatrix(*(src.Entries_[i][j]));
1082 delete Entries_[i][j];
1083 Entries_[i][j] = new Epetra_SerialDenseMatrix(View, Values_ThisBlockEntry,
1084 MyLda, MyRowDim, MyColDim );
1085 }
1086 }
1087 StorageOptimized_ = true ;
1088 }
1089
1090 // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << std::endl ;
1091
1092
1093 /* Work on later...
1094 int i, j;
1095
1096 // The purpose of this routine is to make the block entries in each row contiguous in memory
1097 // so that a single call to GEMV or GEMM call be used to compute an entire block row.
1098
1099
1100 bool Contiguous = true; // Assume contiguous is true
1101 for (i=1; i<NumMyBlockRows_; i++){
1102 int NumEntries = NumBlockEntriesPerRow_[i-1];
1103 int NumAllocatedEntries = NumAllocatedBlockEntriesPerRow_[i-1];
1104
1105 // Check if NumEntries is same as NumAllocatedEntries and
1106 // check if end of beginning of current row starts immediately after end of previous row.
1107 if ((NumEntries!=NumAllocatedEntries) || (Entries_[i]!=Entries_[i-1]+NumEntries)) {
1108 Contiguous = false;
1109 break;
1110 }
1111 }
1112
1113 // NOTE: At the end of the above loop set, there is a possibility that NumEntries and NumAllocatedEntries
1114 // for the last row could be different, but I don't think it matters.
1115
1116
1117 if ((CV_==View) && !Contiguous) EPETRA_CHK_ERR(-1); // This is user data, it's not contiguous and we can't make it so.
1118
1119 int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
1120 if (ierr) EPETRA_CHK_ERR(ierr);
1121
1122 if (Contiguous) return(0); // Everything is done. Return
1123
1124 // Compute Number of Nonzero entries (Done in FillComplete, but we may not have been there yet.)
1125 int numMyNonzeros = Graph_->NumMyNonzeros();
1126
1127 // Allocate one big integer array for all index values
1128 All_Values_ = new double[numMyNonzeros];
1129
1130 // Set Entries_ to point into All_Entries_
1131
1132 double * tmp = All_Values_;
1133 for (i=0; i<NumMyBlockRows_; i++) {
1134 int NumEntries = NumEntriesPerBlockRow_[i];
1135 for (j=0; j<NumEntries; j++) tmp[j] = Entries_[i][j];
1136 if (Entries_[i] !=0) delete [] Entries_[i];
1137 Entries_[i] = tmp;
1138 tmp += NumEntries;
1139 }
1140 */
1141 // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << std::endl ;
1142 return(0);
1143}
1144//==========================================================================
1146 int Length,
1147 int & NumEntries,
1148 double *values,
1149 int * Indices) const
1150{
1151 (void)GlobalRow;
1152 (void)Length;
1153 (void)NumEntries;
1154 (void)values;
1155 (void)Indices;
1156 std::cout << "Must implement..." << std::endl;
1157 return(0);
1158}
1159//==========================================================================
1160int Epetra_VbrMatrix::ExtractGlobalBlockRowPointers(int BlockRow, int maxNumBlockEntries,
1161 int & RowDim, int & NumBlockEntries,
1162 int * BlockIndices,
1163 Epetra_SerialDenseMatrix** & Entries) const {
1164
1165 bool indicesAreLocal = false;
1166 EPETRA_CHK_ERR(ExtractBlockRowPointers(BlockRow, maxNumBlockEntries, RowDim, NumBlockEntries, BlockIndices,
1167 Entries, indicesAreLocal));
1168 return(0);
1169}
1170
1171//==========================================================================
1172int Epetra_VbrMatrix::ExtractMyBlockRowPointers(int BlockRow, int maxNumBlockEntries,
1173 int & RowDim, int & NumBlockEntries,
1174 int * BlockIndices,
1175 Epetra_SerialDenseMatrix** & Entries) const {
1176
1177 bool indicesAreLocal = true;
1178 EPETRA_CHK_ERR(ExtractBlockRowPointers(BlockRow,maxNumBlockEntries , RowDim, NumBlockEntries, BlockIndices,
1179 Entries, indicesAreLocal));
1180 return(0);
1181}
1182
1183//==========================================================================
1184int Epetra_VbrMatrix::ExtractBlockRowPointers(int BlockRow, int maxNumBlockEntries,
1185 int & RowDim, int & NumBlockEntries,
1186 int * BlockIndices,
1187 Epetra_SerialDenseMatrix** & Entries,
1188 bool indicesAreLocal) const {
1189 int ierr = 0;
1190 if (!indicesAreLocal) {
1191 ierr = Graph_->ExtractGlobalRowCopy(BlockRow, maxNumBlockEntries,
1192 NumBlockEntries, BlockIndices);
1193 BlockRow = LRID(BlockRow);
1194 }
1195 else {
1196 ierr = Graph_->ExtractMyRowCopy(BlockRow, maxNumBlockEntries,
1197 NumBlockEntries, BlockIndices);
1198 }
1199 if (ierr) EPETRA_CHK_ERR(ierr);
1200
1201 RowDim = ElementSizeList_[BlockRow];
1202
1203 Entries = Entries_[BlockRow];
1204
1205
1206 EPETRA_CHK_ERR(ierr);
1207 return(0);
1208}
1209
1210//==========================================================================
1211int Epetra_VbrMatrix::BeginExtractGlobalBlockRowCopy(int BlockRow, int maxNumBlockEntries,
1212 int & RowDim, int & NumBlockEntries,
1213 int * BlockIndices, int * ColDims) const {
1214
1215 bool indicesAreLocal = false;
1216 EPETRA_CHK_ERR(BeginExtractBlockRowCopy(BlockRow, maxNumBlockEntries, RowDim, NumBlockEntries, BlockIndices,
1217 ColDims, indicesAreLocal));
1218 return(0);
1219}
1220
1221//==========================================================================
1222int Epetra_VbrMatrix::BeginExtractMyBlockRowCopy(int BlockRow, int maxNumBlockEntries,
1223 int & RowDim, int & NumBlockEntries,
1224 int * BlockIndices, int * ColDims) const {
1225
1226 bool indicesAreLocal = true;
1227 EPETRA_CHK_ERR(BeginExtractBlockRowCopy(BlockRow,maxNumBlockEntries , RowDim, NumBlockEntries, BlockIndices,
1228 ColDims, indicesAreLocal));
1229 return(0);
1230}
1231
1232//==========================================================================
1233int Epetra_VbrMatrix::BeginExtractBlockRowCopy(int BlockRow, int maxNumBlockEntries,
1234 int & RowDim, int & NumBlockEntries,
1235 int * BlockIndices, int * ColDims,
1236 bool indicesAreLocal) const {
1237 int ierr = 0;
1238 if (!indicesAreLocal)
1239 ierr = Graph_->ExtractGlobalRowCopy(BlockRow, maxNumBlockEntries, NumBlockEntries, BlockIndices);
1240 else
1241 ierr = Graph_->ExtractMyRowCopy(BlockRow, maxNumBlockEntries, NumBlockEntries, BlockIndices);
1242 if (ierr) EPETRA_CHK_ERR(ierr);
1243
1244 bool ExtractView = false;
1245 ierr = SetupForExtracts(BlockRow, RowDim, NumBlockEntries, ExtractView, indicesAreLocal);
1246 if (ierr) EPETRA_CHK_ERR(ierr);
1247
1248 EPETRA_CHK_ERR(ExtractBlockDimsCopy(NumBlockEntries, ColDims));
1249 return(0);
1250}
1251
1252//==========================================================================
1253int Epetra_VbrMatrix::SetupForExtracts(int BlockRow, int & RowDim, int NumBlockEntries, bool ExtractView,
1254 bool indicesAreLocal) const {
1255
1256 if (!indicesAreLocal) BlockRow = LRID(BlockRow); // Normalize row range
1257 CurExtractBlockRow_ = BlockRow;
1258 CurExtractEntry_ = 0;
1259 CurExtractNumBlockEntries_ = NumBlockEntries;
1260 CurExtractIndicesAreLocal_ = indicesAreLocal;
1261 CurExtractView_ = ExtractView;
1263 RowDim = CurRowDim_;
1264
1265 return(0);
1266}
1267
1268//==========================================================================
1269int Epetra_VbrMatrix::ExtractBlockDimsCopy(int NumBlockEntries, int * ColDims) const {
1270
1271 for (int i=0; i<NumBlockEntries; i++) {
1272 ColDims[i] = Entries_[CurExtractBlockRow_][i]->N();
1273 }
1274 return(0);
1275}
1276
1277//==========================================================================
1279 double * values,
1280 int LDA,
1281 bool SumInto) const
1282{
1283 (void)SumInto;
1284 if (CurExtractEntry_==-1) EPETRA_CHK_ERR(-1); // No BeginCopy routine was called
1285 int CurColDim = Entries_[CurExtractBlockRow_][CurExtractEntry_]->N();
1286 if (LDA*CurColDim>SizeOfValues) EPETRA_CHK_ERR(-2); // Not enough space
1287
1289 int CurRowDim = CurEntries->M();
1290 int CurLDA = CurEntries->LDA();
1291
1292 CurExtractEntry_++; // Increment Entry Pointer
1293
1294 double* vals = CurEntries->A();
1295 if (LDA==CurRowDim && CurLDA==CurRowDim) // Columns of the entry are contiguous, can treat as single array
1296 for (int ii=0; ii<CurRowDim*CurColDim; ++ii)
1297 values[ii] = vals[ii];
1298 else {
1299 double * CurTargetCol = values;
1300 double * CurSourceCol = vals;
1301
1302 for (int jj=0; jj<CurColDim; jj++) {
1303 for (int ii=0; ii<CurRowDim; ++ii)
1304 CurTargetCol[ii] = CurSourceCol[ii];
1305 CurTargetCol += LDA;
1306 CurSourceCol += CurLDA;
1307 }
1308 }
1309
1310 return(0);
1311}
1312
1313//==========================================================================
1315 int & RowDim, int & NumBlockEntries,
1316 int * & BlockIndices) const
1317{
1318
1319 bool indicesAreLocal = false;
1320 EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries,
1321 BlockIndices, indicesAreLocal));
1322 return(0);
1323}
1324
1325//==========================================================================
1326int Epetra_VbrMatrix::BeginExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
1327 int * & BlockIndices) const
1328{
1329
1330 bool indicesAreLocal = true;
1331 EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries, BlockIndices,
1332 indicesAreLocal));
1333 return(0);
1334}
1335
1336//==========================================================================
1337int Epetra_VbrMatrix::BeginExtractBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
1338 int * & BlockIndices,
1339 bool indicesAreLocal) const
1340{
1341 int ierr = 0;
1342 if (!indicesAreLocal)
1343 ierr = Graph_->ExtractGlobalRowView(BlockRow, NumBlockEntries, BlockIndices);
1344 else
1345 ierr = Graph_->ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices);
1346 if (ierr) EPETRA_CHK_ERR(ierr);
1347
1348 bool ExtractView = true;
1349 ierr = SetupForExtracts(BlockRow, RowDim, NumBlockEntries, ExtractView, indicesAreLocal);
1350 if (ierr) EPETRA_CHK_ERR(ierr);
1351
1352 return(0);
1353}
1354
1355//==========================================================================
1357{
1359
1360 CurExtractEntry_++; // Increment Entry Pointer
1361 return(0);
1362}
1363
1364//==========================================================================
1365int Epetra_VbrMatrix::ExtractGlobalBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
1366 int * & BlockIndices,
1367 Epetra_SerialDenseMatrix** & Entries) const
1368{
1369
1370 Entries = Entries_[LRID(BlockRow)]; // Pointer to Array of pointers for this row's block entries
1371 bool indicesAreLocal = false;
1372 EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries,
1373 BlockIndices,
1374 indicesAreLocal));
1375 return(0);
1376}
1377
1378//==========================================================================
1379int Epetra_VbrMatrix::ExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
1380 int * & BlockIndices,
1381 Epetra_SerialDenseMatrix** & Entries) const
1382{
1383
1384 Entries = Entries_[BlockRow]; // Pointer to Array of pointers for this row's block entries
1385 bool indicesAreLocal = true;
1386 EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries, BlockIndices,
1387 indicesAreLocal));
1388 return(0);
1389}
1390
1391//==============================================================================
1393
1394 if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
1395 if (!RowMap().SameAs(Diagonal.Map())) EPETRA_CHK_ERR(-2); // Maps must be the same
1396 double * diagptr = Diagonal.Values();
1397 for (int i=0; i<NumMyBlockRows_; i++){
1398 int BlockRow = GRID64(i); // FIXME long long
1399 int RowDim = ElementSizeList_[i];
1400 int NumEntries = NumBlockEntriesPerRow_[i];
1401 int * Indices = Indices_[i];
1402 for (int j=0; j<NumEntries; j++) {
1403 int BlockCol = GCID64(Indices[j]);// FIXME long long
1404 if (BlockRow==BlockCol) {
1405 CopyMatDiag(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim,
1406 Entries_[i][j]->N(), diagptr+FirstPointInElementList_[i]);
1407 break;
1408 }
1409 }
1410 }
1411 return(0);
1412}
1413//==============================================================================
1415
1416 if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
1417 if (!RowMap().SameAs(Diagonal.Map())) EPETRA_CHK_ERR(-2); // Maps must be the same
1418 int ierr = 0;
1419 double * diagptr = (double *) Diagonal.Values(); // Dangerous but being lazy
1420 for (int i=0; i<NumMyBlockRows_; i++){
1421 int BlockRow = GRID64(i);// FIXME long long
1422 int RowDim = ElementSizeList_[i];
1423 int NumEntries = NumBlockEntriesPerRow_[i];
1424 int * Indices = Indices_[i];
1425 bool DiagMissing = true;
1426 for (int j=0; j<NumEntries; j++) {
1427 int BlockCol = GCID64(Indices[j]);// FIXME long long
1428 if (BlockRow==BlockCol) {
1429 ReplaceMatDiag(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim, Entries_[i][j]->N(),
1430 diagptr+FirstPointInElementList_[i]);
1431 DiagMissing = false;
1432 break;
1433 }
1434 }
1435 if (DiagMissing) ierr = 1; // flag a warning error
1436 }
1437 NormOne_ = -1.0; // Reset Norm so it will be recomputed.
1438 NormInf_ = -1.0; // Reset Norm so it will be recomputed.
1439 NormFrob_ = -1.0;
1440 EPETRA_CHK_ERR(ierr);
1441 return(0);
1442}
1443//==============================================================================
1444int Epetra_VbrMatrix::BeginExtractBlockDiagonalCopy(int maxNumBlockDiagonalEntries,
1445 int & NumBlockDiagonalEntries, int * RowColDims ) const{
1446
1447 if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
1448 CurBlockDiag_ = 0; // Initialize pointer
1449 NumBlockDiagonalEntries = NumMyBlockRows_;
1450 if (NumBlockDiagonalEntries>maxNumBlockDiagonalEntries) EPETRA_CHK_ERR(-2);
1451 EPETRA_CHK_ERR(RowMap().ElementSizeList(RowColDims));
1452 return(0);
1453}
1454//==============================================================================
1455int Epetra_VbrMatrix::ExtractBlockDiagonalEntryCopy(int SizeOfValues, double * values, int LDA, bool SumInto ) const {
1456
1457 if (CurBlockDiag_==-1) EPETRA_CHK_ERR(-1); // BeginExtractBlockDiagonalCopy was not called
1458 int i = CurBlockDiag_;
1459 int BlockRow = i;
1460 int RowDim = ElementSizeList_[i];
1461 int NumEntries = NumBlockEntriesPerRow_[i];
1462 int * Indices = Indices_[i];
1463 for (int j=0; j<NumEntries; j++) {
1464 int Col = Indices[j];
1465 if (BlockRow==Col) {
1466 int ColDim = Entries_[i][j]->N();
1467 if (LDA*ColDim>SizeOfValues) EPETRA_CHK_ERR(-2); // Not enough room in values
1468 CopyMat(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim, ColDim, values, LDA, SumInto);
1469 break;
1470 }
1471 }
1472 CurBlockDiag_++; // Increment counter
1473 return(0);
1474}
1475//==============================================================================
1476int Epetra_VbrMatrix::BeginExtractBlockDiagonalView(int & NumBlockDiagonalEntries,
1477 int * & RowColDims ) const {
1478
1479 if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
1480 CurBlockDiag_ = 0; // Initialize pointer
1481 NumBlockDiagonalEntries = NumMyBlockRows_;
1482 RowColDims = ElementSizeList_;
1483 return(0);
1484}
1485//==============================================================================
1486int Epetra_VbrMatrix::ExtractBlockDiagonalEntryView(double * & values, int & LDA) const {
1487
1488 if (CurBlockDiag_==-1) EPETRA_CHK_ERR(-1); // BeginExtractBlockDiagonalCopy was not called
1489 int i = CurBlockDiag_;
1490 int BlockRow = i;
1491 int NumEntries = NumBlockEntriesPerRow_[i];
1492 int * Indices = Indices_[i];
1493 for (int j=0; j<NumEntries; j++) {
1494 int Col = Indices[j];
1495 if (BlockRow==Col) {
1496 values = Entries_[i][j]->A();
1497 LDA = Entries_[i][j]->LDA();
1498 break;
1499 }
1500 }
1501 CurBlockDiag_++; // Increment counter
1502 return(0);
1503}
1504//=============================================================================
1505int Epetra_VbrMatrix::CopyMatDiag(double * A, int LDA, int NumRows, int NumCols,
1506 double * Diagonal) const {
1507
1508 int i;
1509 double * ptr1 = Diagonal;
1510 double * ptr2;
1511 int ndiags = EPETRA_MIN(NumRows,NumCols);
1512
1513 for (i=0; i<ndiags; i++) {
1514 ptr2 = A + i*LDA+i;
1515 *ptr1++ = *ptr2;
1516 }
1517 return(0);
1518}
1519//=============================================================================
1520int Epetra_VbrMatrix::ReplaceMatDiag(double * A, int LDA, int NumRows, int NumCols,
1521 double * Diagonal) {
1522
1523 int i;
1524 double * ptr1 = Diagonal;
1525 double * ptr2;
1526 int ndiags = EPETRA_MIN(NumRows,NumCols);
1527
1528 for (i=0; i<ndiags; i++) {
1529 ptr2 = A + i*LDA+i;
1530 *ptr2 = *ptr1++;
1531 }
1532 return(0);
1533}
1534//=============================================================================
1536
1537 int outval = 0;
1538
1539 for (int i=0; i<NumMyBlockRows_; i++){
1540 int NumBlockEntries = NumMyBlockEntries(i);
1541 int NumEntries = 0;
1542 for (int j=0; j<NumBlockEntries; j++) NumEntries += Entries_[i][j]->N();
1543 outval = EPETRA_MAX(outval,NumEntries);
1544 }
1545 return(outval);
1546}
1547//=============================================================================
1548int Epetra_VbrMatrix::NumMyRowEntries(int MyRow, int & NumEntries) const {
1549
1550 int BlockRow, BlockOffset;
1551 int ierr = RowMap().FindLocalElementID(MyRow, BlockRow, BlockOffset); if (ierr!=0) EPETRA_CHK_ERR(ierr);
1552
1553 int NumBlockEntries = NumMyBlockEntries(BlockRow);
1554 NumEntries = 0;
1555 for (int i=0; i<NumBlockEntries; i++) NumEntries += Entries_[BlockRow][i]->N();
1556 return(0);
1557}
1558//=============================================================================
1560 int Length,
1561 int & NumEntries,
1562 double *values,
1563 int * Indices) const
1564{
1565 if (!Filled()) EPETRA_CHK_ERR(-1); // Can't extract row unless matrix is filled
1566 if (!IndicesAreLocal()) EPETRA_CHK_ERR(-2);
1567
1568 int ierr = 0;
1569 int BlockRow, BlockOffset;
1570 ierr = RowMap().FindLocalElementID(MyRow, BlockRow, BlockOffset);
1571 if (ierr!=0) EPETRA_CHK_ERR(ierr);
1572
1573 int RowDim, NumBlockEntries;
1574 int * BlockIndices;
1575 Epetra_SerialDenseMatrix ** ValBlocks;
1576 ierr = ExtractMyBlockRowView(BlockRow, RowDim, NumBlockEntries,
1577 BlockIndices, ValBlocks);
1578 if (ierr!=0) EPETRA_CHK_ERR(ierr);
1579
1580 int * ColFirstPointInElementList = FirstPointInElementList_;
1581 if (Importer()!=0) {
1582 ColFirstPointInElementList = ColMap().FirstPointInElementList();
1583 }
1584 NumEntries = 0;
1585 for (int i=0; i<NumBlockEntries; i++) {
1586 int ColDim = ValBlocks[i]->N();
1587 NumEntries += ColDim;
1588 if (NumEntries>Length) EPETRA_CHK_ERR(-3); // Not enough space
1589 int LDA = ValBlocks[i]->LDA();
1590 double * A = ValBlocks[i]->A()+BlockOffset; // Point to first element in row
1591 int Index = ColFirstPointInElementList[BlockIndices[i]];
1592 for (int j=0; j < ColDim; j++) {
1593 *values++ = *A;
1594 A += LDA;
1595 *Indices++ = Index++;
1596 }
1597 }
1598
1599 return(0);
1600}
1601//=============================================================================
1602int Epetra_VbrMatrix::Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const
1603{
1604 //
1605 // This function forms the product y = A * x or y = A' * x
1606 //
1607
1608 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
1609
1610 int i, j;
1611 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
1612 int * FirstPointInElement = FirstPointInElementList_;
1613 int * ElementSize = ElementSizeList_;
1614 int ** Indices = Indices_;
1616
1617 double * xp = (double*)x.Values();
1618 double *yp = (double*)y.Values();
1619
1620 int * ColElementSizeList = ColMap().ElementSizeList();
1621 int * ColFirstPointInElementList = ColMap().FirstPointInElementList();
1622
1623 bool x_and_y_same = false;
1624 Epetra_Vector* ytemp = 0;
1625 if (xp == yp) {
1626 x_and_y_same = true;
1627 ytemp = new Epetra_Vector(y.Map());
1628 yp = (double*)ytemp->Values();
1629 }
1630
1631 Epetra_MultiVector& yref = x_and_y_same ? *ytemp : y;
1632 //
1633 //Now yref is a reference to ytemp if x_and_y_same == true, otherwise it is a reference
1634 //to y.
1635 //
1636
1637 if (!TransA) {
1638
1639 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
1640 if (Importer()!=0) {
1641 if (ImportVector_!=0) {
1642 if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
1643 }
1644 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
1646 xp = (double*)ImportVector_->Values();
1647 ColElementSizeList = ColMap().ElementSizeList(); // The Import map will always have an existing ElementSizeList
1648 ColFirstPointInElementList = ColMap().FirstPointInElementList(); // Import map will always have an existing ...
1649 }
1650
1651
1652 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
1653 if (Exporter()!=0) {
1654 if (ExportVector_!=0) {
1655 if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
1656 }
1657 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
1658 yp = (double*)ExportVector_->Values();
1659 }
1660
1661 // Do actual computation
1662 int NumMyRows_ = NumMyRows();
1663 for (i=0; i< NumMyRows_; i++) yp[i] = 0.0; // Initialize y
1664
1665 for (i=0; i < NumMyBlockRows_; i++) {
1666 int NumEntries = *NumBlockEntriesPerRow++;
1667 int * BlockRowIndices = *Indices++;
1668 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
1669 double * cury = yp + *FirstPointInElement++;
1670 int RowDim = *ElementSize++;
1671 for (j=0; j < NumEntries; j++) {
1672 //sum += BlockRowValues[j] * xp[BlockRowIndices[j]];
1673 double * A = BlockRowValues[j]->A();
1674 int LDA = BlockRowValues[j]->LDA();
1675 int Index = BlockRowIndices[j];
1676 double * curx = xp + ColFirstPointInElementList[Index];
1677 int ColDim = ColElementSizeList[Index];
1678 GEMV('N', RowDim, ColDim, 1.0, A, LDA, curx, 1.0, cury);
1679 }
1680 }
1681 if (Exporter()!=0) {
1682 yref.PutScalar(0.0);
1683 EPETRA_CHK_ERR(yref.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
1684 }
1685 // Handle case of rangemap being a local replicated map
1686 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(yref.Reduce());
1687 }
1688
1689 else { // Transpose operation
1690
1691
1692 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
1693
1694 if (Exporter()!=0) {
1695 if (ExportVector_!=0) {
1696 if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
1697 }
1698 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
1700 xp = (double*)ExportVector_->Values();
1701 }
1702
1703 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
1704 if (Importer()!=0) {
1705 if (ImportVector_!=0) {
1706 if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
1707 }
1708 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
1709 yp = (double*)ImportVector_->Values();
1710 ColElementSizeList = ColMap().ElementSizeList(); // The Import map will always have an existing ElementSizeList
1711 ColFirstPointInElementList = ColMap().FirstPointInElementList(); // Import map will always have an existing ...
1712 }
1713
1714 // Do actual computation
1715 int NumMyCols_ = NumMyCols();
1716 for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
1717
1718 for (i=0; i < NumMyBlockRows_; i++) {
1719 int NumEntries = *NumBlockEntriesPerRow++;
1720 int * BlockRowIndices = *Indices++;
1721 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
1722 double * curx = xp + *FirstPointInElement++;
1723 int RowDim = *ElementSize++;
1724 for (j=0; j < NumEntries; j++) {
1725 //yp[BlockRowIndices[j]] += BlockRowValues[j] * xp[i];
1726 double * A = BlockRowValues[j]->A();
1727 int LDA = BlockRowValues[j]->LDA();
1728 int Index = BlockRowIndices[j];
1729 double * cury = yp + ColFirstPointInElementList[Index];
1730 int ColDim = ColElementSizeList[Index];
1731 GEMV('T', RowDim, ColDim, 1.0, A, LDA, curx, 1.0, cury);
1732
1733 }
1734 }
1735 if (Importer()!=0) {
1736 yref.PutScalar(0.0); // Make sure target is zero
1737 EPETRA_CHK_ERR(yref.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
1738 }
1739 // Handle case of rangemap being a local replicated map
1740 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(yref.Reduce());
1741 }
1742
1743 if (x_and_y_same == true) {
1744 y = *ytemp;
1745 delete ytemp;
1746 }
1747
1749 return(0);
1750}
1751
1752//=============================================================================
1754 //
1755 // This function forms the product Y = A * Y or Y = A' * X
1756 //
1757
1758 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
1759
1760 int i;
1761 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
1762 int ** Indices = Indices_;
1764
1765 int * RowElementSizeList = ElementSizeList_;
1766 int * RowFirstPointInElementList = FirstPointInElementList_;
1767 int * ColElementSizeList = ColMap().ElementSizeList();
1768 int * ColFirstPointInElementList = ColMap().FirstPointInElementList();
1769
1770
1771 int NumVectors = X.NumVectors();
1772 double **Xp = (double**)X.Pointers();
1773 double **Yp = (double**)Y.Pointers();
1774
1775 bool x_and_y_same = false;
1776 Epetra_MultiVector* Ytemp = 0;
1777 if (*Xp == *Yp) {
1778 x_and_y_same = true;
1779 Ytemp = new Epetra_MultiVector(Y.Map(), NumVectors);
1780 Yp = (double**)Ytemp->Pointers();
1781 }
1782
1783 Epetra_MultiVector& Yref = x_and_y_same ? *Ytemp : Y;
1784 //
1785 //Now Yref is a reference to Ytemp if x_and_y_same == true, otherwise it is a reference
1786 //to Y.
1787 //
1788
1789 if (!TransA) {
1790
1791 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
1792 if (Importer()!=0) {
1793 if (ImportVector_!=0) {
1794 if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
1795 }
1796 // Create import vector if needed
1797 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors);
1798
1800 Xp = (double**)ImportVector_->Pointers();
1801 ColElementSizeList = ColMap().ElementSizeList();
1802 ColFirstPointInElementList = ColMap().FirstPointInElementList();
1803 }
1804
1805 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
1806 if (Exporter()!=0) {
1807 if (ExportVector_!=0) {
1808 if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
1809 }
1810 // Create Export vector if needed
1811 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors);
1812
1813 ExportVector_->PutScalar(0.0); // Zero y values
1814 Yp = (double**)ExportVector_->Pointers();
1815 RowElementSizeList = ColMap().ElementSizeList();
1816 RowFirstPointInElementList = ColMap().FirstPointInElementList();
1817 }
1818 else {
1819 Yref.PutScalar(0.0); // Zero y values
1820 }
1821
1822 // Do actual computation
1823 if ( All_Values_ != 0 ) { // Contiguous Storage
1824 if ( ! TransA && *RowElementSizeList <= 4 ) {
1825
1826 int RowDim = *RowElementSizeList ;
1827
1828 Epetra_SerialDenseMatrix* Asub = **Entries;
1829 double *A = Asub->A_ ;
1830
1831 if ( NumVectors == 1 ) {
1832
1833 for (i=0; i < NumMyBlockRows_; i++) {
1834 int NumEntries = *NumBlockEntriesPerRow++;
1835 int * BlockRowIndices = *Indices++;
1836
1837 double * xptr = Xp[0];
1838 double y0 = 0.0;
1839 double y1 = 0.0;
1840 double y2 = 0.0;
1841 double y3 = 0.0;
1842 switch(RowDim) {
1843 case 1:
1844 for (int j=0; j < NumEntries; ++j) {
1845 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1846
1847 y0 += A[0]*xptr[xoff];
1848 A += 1;
1849 }
1850 break;
1851
1852 case 2:
1853 for (int j=0; j < NumEntries; ++j) {
1854 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1855 y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
1856 y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
1857 A += 4 ;
1858 }
1859 break;
1860
1861 case 3:
1862 for (int j=0; j < NumEntries; ++j) {
1863 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1864 y0 += A[0]*xptr[xoff+0] + A[3]*xptr[xoff+1] + A[6]*xptr[xoff+2];
1865 y1 += A[1]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[7]*xptr[xoff+2];
1866 y2 += A[2]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[8]*xptr[xoff+2];
1867 A += 9 ;
1868 }
1869 break;
1870
1871 case 4:
1872 for (int j=0; j < NumEntries; ++j) {
1873 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1874 y0 += A[0]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[8]*xptr[xoff+2] +A[12]*xptr[xoff+3];
1875 y1 += A[1]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[9]*xptr[xoff+2] +A[13]*xptr[xoff+3];
1876 y2 += A[2]*xptr[xoff+0] + A[6]*xptr[xoff+1] + A[10]*xptr[xoff+2] +A[14]*xptr[xoff+3];
1877 y3 += A[3]*xptr[xoff+0] + A[7]*xptr[xoff+1] + A[11]*xptr[xoff+2] +A[15]*xptr[xoff+3];
1878 A += 16 ;
1879 }
1880 break;
1881 default:
1882 assert(false);
1883 }
1884 int yoff = *RowFirstPointInElementList++;
1885 switch(RowDim) {
1886 case 4:
1887 *(Yp[0]+yoff+3) = y3;
1888 case 3:
1889 *(Yp[0]+yoff+2) = y2;
1890 case 2:
1891 *(Yp[0]+yoff+1) = y1;
1892 case 1:
1893 *(Yp[0]+yoff) = y0;
1894 }
1895 }
1896 } else { // NumVectors != 1
1897 double *OrigA = A ;
1898 for (i=0; i < NumMyBlockRows_; i++) {
1899 int NumEntries = *NumBlockEntriesPerRow++;
1900 int yoff = *RowFirstPointInElementList++;
1901 int * BRI = *Indices++;
1902
1903 for (int k=0; k<NumVectors; k++) {
1904 int * BlockRowIndices = BRI;
1905 A = OrigA ;
1906 double * xptr = Xp[k];
1907 double y0 = 0.0;
1908 double y1 = 0.0;
1909 double y2 = 0.0;
1910 double y3 = 0.0;
1911 switch(RowDim) {
1912 case 1:
1913 for (int j=0; j < NumEntries; ++j) {
1914 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1915
1916 y0 += A[0]*xptr[xoff];
1917 A += 1;
1918 }
1919 break;
1920
1921 case 2:
1922 for (int j=0; j < NumEntries; ++j) {
1923 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1924 y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
1925 y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
1926 A += 4 ;
1927 }
1928 break;
1929
1930 case 3:
1931 for (int j=0; j < NumEntries; ++j) {
1932 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1933 y0 += A[0]*xptr[xoff+0] + A[3]*xptr[xoff+1] + A[6]*xptr[xoff+2];
1934 y1 += A[1]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[7]*xptr[xoff+2];
1935 y2 += A[2]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[8]*xptr[xoff+2];
1936 A += 9 ;
1937 }
1938 break;
1939 case 4:
1940 for (int j=0; j < NumEntries; ++j) {
1941 int xoff = ColFirstPointInElementList[*BlockRowIndices++];
1942 y0 += A[0]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[8]*xptr[xoff+2] +A[12]*xptr[xoff+3];
1943 y1 += A[1]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[9]*xptr[xoff+2] +A[13]*xptr[xoff+3];
1944 y2 += A[2]*xptr[xoff+0] + A[6]*xptr[xoff+1] + A[10]*xptr[xoff+2] +A[14]*xptr[xoff+3];
1945 y3 += A[3]*xptr[xoff+0] + A[7]*xptr[xoff+1] + A[11]*xptr[xoff+2] +A[15]*xptr[xoff+3];
1946 A += 16 ;
1947 }
1948 break;
1949 default:
1950 assert(false);
1951 }
1952 switch(RowDim) {
1953 case 4:
1954 *(Yp[k]+yoff+3) = y3;
1955 case 3:
1956 *(Yp[k]+yoff+2) = y2;
1957 case 2:
1958 *(Yp[k]+yoff+1) = y1;
1959 case 1:
1960 *(Yp[k]+yoff) = y0;
1961 }
1962 }
1963 OrigA = A ;
1964 }
1965 } // end else NumVectors != 1
1966
1967 } else {
1968
1969 for (i=0; i < NumMyBlockRows_; i++) {
1970 int NumEntries = *NumBlockEntriesPerRow++;
1971 int * BlockRowIndices = *Indices++;
1972 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
1973 int yoff = *RowFirstPointInElementList++;
1974 int RowDim = *RowElementSizeList++;
1975
1976 FastBlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff,
1977 ColFirstPointInElementList, ColElementSizeList,
1978 BlockRowValues, Xp, Yp, NumVectors);
1979 }
1980 }
1981 } else {
1982 for (i=0; i < NumMyBlockRows_; i++) {
1983 int NumEntries = *NumBlockEntriesPerRow++;
1984 int * BlockRowIndices = *Indices++;
1985 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
1986 int yoff = *RowFirstPointInElementList++;
1987 int RowDim = *RowElementSizeList++;
1988 BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff,
1989 ColFirstPointInElementList, ColElementSizeList,
1990 BlockRowValues, Xp, Yp, NumVectors);
1991 }
1992 }
1993
1994 if (Exporter()!=0) {
1995 Yref.PutScalar(0.0);
1996 EPETRA_CHK_ERR(Yref.Export(*ExportVector_, *Exporter(), Add)); // Fill Y with Values from export vector
1997 }
1998 // Handle case of rangemap being a local replicated map
1999 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Yref.Reduce());
2000 }
2001 else { // Transpose operation
2002
2003
2004 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
2005
2006 if (Exporter()!=0) {
2007 if (ExportVector_!=0) {
2008 if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
2009 }
2010 // Create Export vector if needed
2011 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors);
2012
2014 Xp = (double**)ExportVector_->Pointers();
2015 ColElementSizeList = RowMap().ElementSizeList();
2016 ColFirstPointInElementList = RowMap().FirstPointInElementList();
2017 }
2018
2019 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
2020 if (Importer()!=0) {
2021 if (ImportVector_!=0) {
2022 if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
2023 }
2024 // Create import vector if needed
2025 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors);
2026
2027 ImportVector_->PutScalar(0.0); // Zero y values
2028 Yp = (double**)ImportVector_->Pointers();
2029 RowElementSizeList = ColMap().ElementSizeList();
2030 RowFirstPointInElementList = ColMap().FirstPointInElementList();
2031 }
2032 else
2033 Yref.PutScalar(0.0); // Zero y values
2034
2035 // Do actual computation
2036
2037 if ( All_Values_ != 0 ) { // Contiguous Storage
2038 for (i=0; i < NumMyBlockRows_; i++) {
2039 int NumEntries = *NumBlockEntriesPerRow++;
2040 int * BlockRowIndices = *Indices++;
2041 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2042 int xoff = *RowFirstPointInElementList++;
2043 int RowDim = *RowElementSizeList++;
2044 FastBlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, xoff,
2045 ColFirstPointInElementList, ColElementSizeList,
2046 BlockRowValues, Xp, Yp, NumVectors);
2047 }
2048 } else {
2049 for (i=0; i < NumMyBlockRows_; i++) {
2050 int NumEntries = *NumBlockEntriesPerRow++;
2051 int * BlockRowIndices = *Indices++;
2052 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2053 int xoff = *RowFirstPointInElementList++;
2054 int RowDim = *RowElementSizeList++;
2055 BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, xoff,
2056 ColFirstPointInElementList, ColElementSizeList,
2057 BlockRowValues, Xp, Yp, NumVectors);
2058 }
2059 }
2060
2061 if (Importer()!=0) {
2062 Yref.PutScalar(0.0); // Make sure target is zero
2063 EPETRA_CHK_ERR(Yref.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
2064 }
2065 // Handle case of rangemap being a local replicated map
2066 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Yref.Reduce());
2067 }
2068
2069 UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
2070
2071 if (x_and_y_same == true) {
2072 Y = *Ytemp;
2073 delete Ytemp;
2074 }
2075
2076 return(0);
2077}
2078//=============================================================================
2080 int RowDim,
2081 int NumEntries,
2082 int * BlockIndices,
2083 int RowOff,
2084 int * FirstPointInElementList,
2085 int * ElementSizeList,
2086 double Alpha,
2088 double ** X,
2089 double Beta,
2090 double ** Y,
2091 int NumVectors) const
2092{
2093 //This overloading of BlockRowMultiply is the same as the one below, except
2094 //that this one accepts Alpha and Beta arguments. This BlockRowMultiply is
2095 //called from within the 'solve' methods.
2096 int j, k;
2097 if (!TransA) {
2098 for (j=0; j < NumEntries; j++) {
2099 Epetra_SerialDenseMatrix* Asub = As[j];
2100 double * A = Asub->A();
2101 int LDA = Asub->LDA();
2102 int BlockIndex = BlockIndices[j];
2103 int xoff = FirstPointInElementList[BlockIndex];
2104 int ColDim = ElementSizeList[BlockIndex];
2105
2106 for (k=0; k<NumVectors; k++) {
2107 double * curx = X[k] + xoff;
2108 double * cury = Y[k] + RowOff;
2109
2110 GEMV('N', RowDim, ColDim, Alpha, A, LDA, curx, Beta, cury);
2111 }//end for (k
2112 }//end for (j
2113 }
2114 else { //TransA == true
2115 for (j=0; j < NumEntries; j++) {
2116 double * A = As[j]->A();
2117 int LDA = As[j]->LDA();
2118 int BlockIndex = BlockIndices[j];
2119 int yoff = FirstPointInElementList[BlockIndex];
2120 int ColDim = ElementSizeList[BlockIndex];
2121 for (k=0; k<NumVectors; k++) {
2122 double * curx = X[k] + RowOff;
2123 double * cury = Y[k] + yoff;
2124 GEMV('T', RowDim, ColDim, Alpha, A, LDA, curx, Beta, cury);
2125 }
2126 }
2127 }
2128
2129 return;
2130}
2131//=============================================================================
2133 int RowDim,
2134 int NumEntries,
2135 int * BlockRowIndices,
2136 int yoff,
2137 int * ColFirstPointInElementList,
2138 int * ColElementSizeList,
2139 Epetra_SerialDenseMatrix** BlockRowValues,
2140 double ** Xp,
2141 double ** Yp,
2142 int NumVectors) const
2143{
2144 int j, k;
2145 if (!TransA) {
2146 for (k=0; k<NumVectors; k++) {
2147 double * y = Yp[k] + yoff;
2148 double * xptr = Xp[k];
2149
2150 Epetra_SerialDenseMatrix* Asub = BlockRowValues[0];
2151 double *OrigA = Asub->A_;
2152 int LDA = Asub->LDA_;
2153 int OrigBlockIndex = BlockRowIndices[0];
2154 int ColDim = ColElementSizeList[OrigBlockIndex];
2155
2156 assert( RowDim == ColDim ) ;
2157 assert( RowDim == LDA ) ;
2158
2159 switch(RowDim) {
2160#if 0
2161 case 1:
2162 double y0 = y[0];
2163 double y1 = y[0];
2164 double *A = OrigA ;
2165 for (j=0; j < NumEntries; ++j) {
2166
2167 int BlockIndex = BlockRowIndices[j];
2168 int xoff = ColFirstPointInElementList[BlockIndex];
2169
2170 double *A = OrigA + j * ColDim * LDA ;
2171 double *x = xptr + xoff;
2172 y[0] += A[0]*x[0];
2173 }//end for (j
2174 break;
2175
2176 case 2:
2177 double y0 = y[0];
2178 double y1 = y[0];
2179 double *A = OrigA ;
2180 for (j=0; j < NumEntries; ++j) {
2181
2182 int BlockIndex = BlockRowIndices[j];
2183 int xoff = ColFirstPointInElementList[BlockIndex];
2184
2185 y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
2186 y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
2187 A += ColDim * LDA ;
2188 }
2189 y[0] = y0;
2190 y[1] = y1;
2191 break;
2192
2193 case 3:
2194 for (j=0; j < NumEntries; ++j) {
2195
2196 int BlockIndex = BlockRowIndices[j];
2197 int xoff = ColFirstPointInElementList[BlockIndex];
2198
2199 double *A = OrigA + j * ColDim * LDA ;
2200 double *x = xptr + xoff;
2201 y[0] += A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
2202 y[1] += A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
2203 y[2] += A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
2204 }
2205 break;
2206
2207 case 4:
2208 for (j=0; j < NumEntries; ++j) {
2209
2210 int BlockIndex = BlockRowIndices[j];
2211 int xoff = ColFirstPointInElementList[BlockIndex];
2212
2213 double *A = OrigA + j * ColDim * LDA ;
2214 double *x = xptr + xoff;
2215 y[0] += A[0]*x[0] + A[4]*x[1] + A[8]*x[2] + A[12]*x[3];
2216 y[1] += A[1]*x[0] + A[5]*x[1] + A[9]*x[2] + A[13]*x[3];
2217 y[2] += A[2]*x[0] + A[6]*x[1] + A[10]*x[2] + A[14]*x[3];
2218 y[3] += A[3]*x[0] + A[7]*x[1] + A[11]*x[2] + A[15]*x[3];
2219 }
2220 break;
2221#endif
2222 case 5:
2223 for (j=0; j < NumEntries; ++j) {
2224
2225 int BlockIndex = BlockRowIndices[j];
2226 int xoff = ColFirstPointInElementList[BlockIndex];
2227
2228 double *A = OrigA + j * ColDim * LDA ;
2229 double *x = xptr + xoff;
2230 y[0] += A[0]*x[0] + A[5]*x[1] + A[10]*x[2] + A[15]*x[3] + A[20]*x[4];
2231 y[1] += A[1]*x[0] + A[6]*x[1] + A[11]*x[2] + A[16]*x[3] + A[21]*x[4];
2232 y[2] += A[2]*x[0] + A[7]*x[1] + A[12]*x[2] + A[17]*x[3] + A[22]*x[4];
2233 y[3] += A[3]*x[0] + A[8]*x[1] + A[13]*x[2] + A[18]*x[3] + A[23]*x[4];
2234 y[4] += A[4]*x[0] + A[9]*x[1] + A[14]*x[2] + A[19]*x[3] + A[24]*x[4];
2235 }
2236 break;
2237
2238 case 6:
2239 for (j=0; j < NumEntries; ++j) {
2240
2241 int BlockIndex = BlockRowIndices[j];
2242 int xoff = ColFirstPointInElementList[BlockIndex];
2243
2244 double *A = OrigA + j * ColDim * LDA ;
2245 double *x = xptr + xoff;
2246 y[0] += A[0]*x[0] + A[6]*x[1] + A[12]*x[2] + A[18]*x[3] + A[24]*x[4]
2247 + A[30]*x[5];
2248 y[1] += A[1]*x[0] + A[7]*x[1] + A[13]*x[2] + A[19]*x[3] + A[25]*x[4]
2249 + A[31]*x[5];
2250 y[2] += A[2]*x[0] + A[8]*x[1] + A[14]*x[2] + A[20]*x[3] + A[26]*x[4]
2251 + A[32]*x[5];
2252 y[3] += A[3]*x[0] + A[9]*x[1] + A[15]*x[2] + A[21]*x[3] + A[27]*x[4]
2253 + A[33]*x[5];
2254 y[4] += A[4]*x[0] + A[10]*x[1] + A[16]*x[2] + A[22]*x[3] + A[28]*x[4]
2255 + A[34]*x[5];
2256 y[5] += A[5]*x[0] + A[11]*x[1] + A[17]*x[2] + A[23]*x[3] + A[29]*x[4]
2257 + A[35]*x[5];
2258 }
2259 break;
2260
2261 default:
2262 for (j=0; j < NumEntries; ++j) {
2263
2264 int BlockIndex = BlockRowIndices[j];
2265 int xoff = ColFirstPointInElementList[BlockIndex];
2266
2267 double *A = OrigA + j * ColDim * LDA ;
2268 double *x = xptr + xoff;
2269 GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
2270 }
2271 }//end switch
2272
2273 }//end for (k
2274 }
2275 else { //TransA == true
2276 for (j=0; j < NumEntries; j++) {
2277 double * A = BlockRowValues[j]->A();
2278 int LDA = BlockRowValues[j]->LDA();
2279 int BlockIndex = BlockRowIndices[j];
2280 int Yyoff = ColFirstPointInElementList[BlockIndex];
2281 int ColDim = ColElementSizeList[BlockIndex];
2282 for (k=0; k<NumVectors; k++) {
2283 double * x = Xp[k] + yoff;
2284 double * y = Yp[k] + Yyoff;
2285 GEMV('T', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
2286 }
2287 }
2288 }
2289
2290}
2292 int RowDim,
2293 int NumEntries,
2294 int * BlockIndices,
2295 int RowOff,
2296 int * FirstPointInElementList,
2297 int * ElementSizeList,
2299 double ** X,
2300 double ** Y,
2301 int NumVectors) const
2302{
2303 //This overloading of BlockRowMultiply is the same as the one above, except
2304 //that this one doesn't accept the Alpha and Beta arguments (it assumes that
2305 //they are both 1.0) and contains some inlined unrolled code for certain
2306 //cases (certain block-sizes) rather than calling GEMV every time. This
2307 //BlockRowMultiply is called from within the 'Multiply' methods.
2308 //Note: Scott Hutchinson's Aztec method 'dvbr_sparax_basic' was consulted in
2309 //the optimizing of this method.
2310
2311 int j, k;
2312 if (!TransA) {
2313 for (k=0; k<NumVectors; k++) {
2314 double * y = Y[k] + RowOff;
2315 double * xptr = X[k];
2316
2317 for (j=0; j < NumEntries; ++j) {
2318 Epetra_SerialDenseMatrix* Asub = As[j];
2319 double * A = Asub->A_;
2320 int LDA = Asub->LDA_;
2321 int BlockIndex = BlockIndices[j];
2322 int xoff = FirstPointInElementList[BlockIndex];
2323 int ColDim = ElementSizeList[BlockIndex];
2324
2325 double * x = xptr + xoff;
2326
2327 //Call GEMV if sub-block is non-square or if LDA != RowDim.
2328 if (LDA != RowDim || ColDim != RowDim) {
2329 GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
2330 continue;
2331 }
2332
2333 //It is a big performance win to use inlined, unrolled code for small
2334 //common block sizes rather than calling GEMV.
2335
2336 switch(RowDim) {
2337 case 1:
2338 y[0] += A[0]*x[0];
2339 break;
2340
2341 case 2:
2342 y[0] += A[0]*x[0] + A[2]*x[1];
2343 y[1] += A[1]*x[0] + A[3]*x[1];
2344 break;
2345
2346 case 3:
2347 y[0] += A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
2348 y[1] += A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
2349 y[2] += A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
2350 break;
2351
2352 case 4:
2353 y[0] += A[0]*x[0] + A[4]*x[1] + A[8]*x[2] + A[12]*x[3];
2354 y[1] += A[1]*x[0] + A[5]*x[1] + A[9]*x[2] + A[13]*x[3];
2355 y[2] += A[2]*x[0] + A[6]*x[1] + A[10]*x[2] + A[14]*x[3];
2356 y[3] += A[3]*x[0] + A[7]*x[1] + A[11]*x[2] + A[15]*x[3];
2357 break;
2358
2359 case 5:
2360 y[0] += A[0]*x[0] + A[5]*x[1] + A[10]*x[2] + A[15]*x[3] + A[20]*x[4];
2361 y[1] += A[1]*x[0] + A[6]*x[1] + A[11]*x[2] + A[16]*x[3] + A[21]*x[4];
2362 y[2] += A[2]*x[0] + A[7]*x[1] + A[12]*x[2] + A[17]*x[3] + A[22]*x[4];
2363 y[3] += A[3]*x[0] + A[8]*x[1] + A[13]*x[2] + A[18]*x[3] + A[23]*x[4];
2364 y[4] += A[4]*x[0] + A[9]*x[1] + A[14]*x[2] + A[19]*x[3] + A[24]*x[4];
2365 break;
2366
2367 case 6:
2368 y[0] += A[0]*x[0] + A[6]*x[1] + A[12]*x[2] + A[18]*x[3] + A[24]*x[4]
2369 + A[30]*x[5];
2370 y[1] += A[1]*x[0] + A[7]*x[1] + A[13]*x[2] + A[19]*x[3] + A[25]*x[4]
2371 + A[31]*x[5];
2372 y[2] += A[2]*x[0] + A[8]*x[1] + A[14]*x[2] + A[20]*x[3] + A[26]*x[4]
2373 + A[32]*x[5];
2374 y[3] += A[3]*x[0] + A[9]*x[1] + A[15]*x[2] + A[21]*x[3] + A[27]*x[4]
2375 + A[33]*x[5];
2376 y[4] += A[4]*x[0] + A[10]*x[1] + A[16]*x[2] + A[22]*x[3] + A[28]*x[4]
2377 + A[34]*x[5];
2378 y[5] += A[5]*x[0] + A[11]*x[1] + A[17]*x[2] + A[23]*x[3] + A[29]*x[4]
2379 + A[35]*x[5];
2380 break;
2381
2382 default:
2383 GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
2384 }//end switch
2385 }//end for (k
2386 }//end for (j
2387 }
2388 else { //TransA == true
2389 for (j=0; j < NumEntries; j++) {
2390 double * A = As[j]->A();
2391 int LDA = As[j]->LDA();
2392 int BlockIndex = BlockIndices[j];
2393 int yoff = FirstPointInElementList[BlockIndex];
2394 int ColDim = ElementSizeList[BlockIndex];
2395 for (k=0; k<NumVectors; k++) {
2396 double * x = X[k] + RowOff;
2397 double * y = Y[k] + yoff;
2398 GEMV('T', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
2399 }
2400 }
2401 }
2402
2403 return;
2404}
2405//=============================================================================
2406int Epetra_VbrMatrix::DoSolve(bool Upper, bool TransA,
2407 bool UnitDiagonal,
2408 const Epetra_MultiVector& X,
2409 Epetra_MultiVector& Y) const
2410{
2411 (void)UnitDiagonal;
2412 //
2413 // This function find Y such that LY = X or UY = X or the transpose cases.
2414 //
2415
2416 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
2417
2418 if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
2419 if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
2420 if (!NoDiagonal()) EPETRA_CHK_ERR(-4); // We must use UnitDiagonal
2421
2422 int i;
2423 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2424 int * FirstPointInElement = FirstPointInElementList_;
2425 int * ElementSize = ElementSizeList_;
2426 int ** Indices = Indices_;
2428
2429 int * ColElementSizeList = ElementSizeList_;
2430 int * ColFirstPointInElementList = FirstPointInElementList_;
2431
2432 // If upper, point to last row
2433 if (Upper) {
2434 NumBlockEntriesPerRow += NumMyBlockRows_-1;
2435 FirstPointInElement += NumMyBlockRows_-1;
2436 ElementSize += NumMyBlockRows_-1;
2437 Indices += NumMyBlockRows_-1;
2438 Entries += NumMyBlockRows_-1;
2439 }
2440
2441 double **Yp = (double**)Y.Pointers();
2442
2443 int NumVectors = X.NumVectors();
2444
2445 if (X.Pointers() != Y.Pointers()) Y = X; // Copy X into Y if they are not the same multivector
2446
2447 bool Case1 = (((!TransA) && Upper) || (TransA && !Upper));
2448 // Case 2 = (((TransA) && Upper) || (!TransA) && Lower);
2449 if (Case1) {
2450 for (i=0; i < NumMyBlockRows_; i++) {
2451 int NumEntries = *NumBlockEntriesPerRow--;
2452 int * BlockRowIndices = *Indices--;
2453 Epetra_SerialDenseMatrix** BlockRowValues = *Entries--;
2454 int yoff = *FirstPointInElement--;
2455 int RowDim = *ElementSize--;
2456 BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff,
2457 ColFirstPointInElementList, ColElementSizeList,
2458 1.0, BlockRowValues, Yp, -1.0, Yp, NumVectors);
2459 }
2460 }
2461 else {
2462 for (i=0; i < NumMyBlockRows_; i++) {
2463 int NumEntries = *NumBlockEntriesPerRow++;
2464 int * BlockRowIndices = *Indices++;
2465 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2466 int yoff = *FirstPointInElement++;
2467 int RowDim = *ElementSize++;
2468 BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff,
2469 ColFirstPointInElementList, ColElementSizeList,
2470 1.0, BlockRowValues, Yp, -1.0, Yp, NumVectors);
2471 }
2472 }
2473
2474 UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
2475 return(0);
2476}
2477
2478//=============================================================================
2480 //
2481 // Put inverse of the sum of absolute values of the ith row of A in x[i].
2482 //
2483 EPETRA_CHK_ERR(InverseSums(true, x));
2484 return(0);
2485}
2486
2487//=============================================================================
2489 //
2490 // Put inverse of the sum of absolute values of the jth column of A in x[j].
2491 //
2492 EPETRA_CHK_ERR(InverseSums(false, x));
2493 return(0);
2494}
2495//=============================================================================
2497 //
2498 // Put inverse of the sum of absolute values of the ith row of A in x[i].
2499 //
2500
2501 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
2502 bool hasOperatorMap = false;
2503 if (DoRows) {
2504 if ( !Graph().RangeMap().SameAs(x.Map()) ) {
2505 hasOperatorMap = OperatorRangeMap().SameAs(x.Map());
2506 if( !hasOperatorMap )
2507 EPETRA_CHK_ERR(-2);
2508 }
2509 }
2510 else {
2511 if ( !Graph().DomainMap().SameAs(x.Map()) ) {
2512 hasOperatorMap = OperatorDomainMap().SameAs(x.Map());
2513 if( !hasOperatorMap )
2514 EPETRA_CHK_ERR(-2);
2515 }
2516 }
2517 int ierr = 0;
2518 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2519 int ** Indices = Indices_;
2521
2522 int * RowElementSizeList = ElementSizeList_;
2523 int * RowFirstPointInElementList = FirstPointInElementList_;
2524 int * ColElementSizeList = ElementSizeList_;
2525 int * ColFirstPointInElementList = FirstPointInElementList_;
2526 if (Importer()!=0) {
2527 ColElementSizeList = ColMap().ElementSizeList();
2528 ColFirstPointInElementList = ColMap().FirstPointInElementList();
2529 }
2530
2531 x.PutScalar(0.0); // Zero out result vector
2532
2533 double * xp = (double*)x.Values();
2534
2535 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
2536 Epetra_Vector * x_tmp = 0;
2537 if (!DoRows) {
2538 if (Importer()!=0) {
2539 x_tmp = new Epetra_Vector(ColMap()); // Create import vector if needed
2540 xp = (double*)x_tmp->Values();
2541 }
2542 }
2543
2544 for (int i=0; i < NumMyBlockRows_; i++) {
2545 int NumEntries = *NumBlockEntriesPerRow++;
2546 int * BlockRowIndices = *Indices++;
2547 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2548 int xoff = *RowFirstPointInElementList++;
2549 int RowDim = *RowElementSizeList++;
2550 if (DoRows) {
2551 for (int ii=0; ii < NumEntries; ii++) {
2552 double * xptr = xp+xoff;
2553 double * A = BlockRowValues[ii]->A();
2554 int LDA = BlockRowValues[ii]->LDA();
2555 int BlockIndex = BlockRowIndices[ii];
2556 int ColDim = ColElementSizeList[BlockIndex];
2557 for (int j=0; j<ColDim; j++) {
2558 double * curEntry = A + j*LDA;
2559 for (int k=0; k<RowDim; k++)
2560 xptr[k] += std::abs(*curEntry++);
2561 }
2562 }
2563 }
2564 else {
2565 for (int ii=0; ii < NumEntries; ii++) {
2566 double * A = BlockRowValues[ii]->A();
2567 int LDA = BlockRowValues[ii]->LDA();
2568 int BlockIndex = BlockRowIndices[ii];
2569 int off = ColFirstPointInElementList[BlockIndex];
2570 int ColDim = ColElementSizeList[BlockIndex];
2571 double * curx = xp+off;
2572 for (int j=0; j<ColDim; j++) {
2573 double * curEntry = A + j*LDA;
2574 for (int k=0; k<RowDim; k++)
2575 curx[j] += std::abs(*curEntry++);
2576 }
2577 }
2578 }
2579 }
2580
2581 if (!DoRows) {
2582 if (Importer()!=0){
2583 Epetra_Vector *x_blocked = 0;
2584 if(hasOperatorMap)
2585 x_blocked = new Epetra_Vector( ::View, DomainMap(), &x[0] );
2586 else
2587 x_blocked = &x;
2588 x_blocked->PutScalar(0.0);
2589 EPETRA_CHK_ERR(x_blocked->Export(*x_tmp, *Importer(), Add)); // Fill x with Values from import vector
2590 if(hasOperatorMap)
2591 delete x_blocked;
2592 delete x_tmp;
2593 xp = (double*) x.Values();
2594 }
2595 }
2596 int NumMyRows_ = NumMyRows();
2597 {for (int i=0; i < NumMyRows_; i++) {
2598 double scale = xp[i];
2599 if (scale<Epetra_MinDouble) {
2600 if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero row/col sum found (supercedes ierr = 2)
2601 else if (ierr!=1) ierr = 2;
2602 xp[i] = Epetra_MaxDouble;
2603 }
2604 else
2605 xp[i] = 1.0/scale;
2606 }}
2608
2609 EPETRA_CHK_ERR(ierr);
2610 return(0);
2611}
2612//=============================================================================
2614 //
2615 // Multiply the ith row of A by x[i].
2616 //
2617 EPETRA_CHK_ERR(Scale(true, x));
2618 return(0);
2619}
2620
2621//=============================================================================
2623 //
2624 // Multiply the jth column of A by x[j].
2625 //
2626 EPETRA_CHK_ERR(Scale (false, x));
2627 return(0);
2628}
2629//=============================================================================
2630int Epetra_VbrMatrix::Scale(bool DoRows, const Epetra_Vector& x) {
2631
2632 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
2633 bool hasOperatorMap = false;
2634 if (DoRows) {
2635 if ( !Graph().RangeMap().SameAs(x.Map()) ) {
2636 hasOperatorMap = OperatorRangeMap().SameAs(x.Map());
2637 if( !hasOperatorMap )
2638 EPETRA_CHK_ERR(-2);
2639 }
2640 }
2641 else {
2642 if ( !Graph().DomainMap().SameAs(x.Map()) ) {
2643 hasOperatorMap = OperatorDomainMap().SameAs(x.Map());
2644 if( !hasOperatorMap )
2645 EPETRA_CHK_ERR(-2);
2646 }
2647 }
2648 int ierr = 0;
2649 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2650 int ** Indices = Indices_;
2652
2653 int * RowElementSizeList = ElementSizeList_;
2654 int * RowFirstPointInElementList = FirstPointInElementList_;
2655 int * ColElementSizeList = ElementSizeList_;
2656 int * ColFirstPointInElementList = FirstPointInElementList_;
2657 if (Importer()!=0) {
2658 ColElementSizeList = ColMap().ElementSizeList();
2659 ColFirstPointInElementList = ColMap().FirstPointInElementList();
2660 }
2661
2662 double * xp = (double*)x.Values();
2663
2664 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
2665 Epetra_Vector * x_tmp = 0;
2666 if (!DoRows) {
2667 if (Importer()!=0) {
2668 Epetra_Vector *x_blocked = 0;
2669 if(hasOperatorMap)
2670 x_blocked = new Epetra_Vector( ::View, Graph().DomainMap(), (double *) &x[0] );
2671 else
2672 x_blocked = (Epetra_Vector *) &x;
2673 x_tmp = new Epetra_Vector(ColMap()); // Create import vector if needed
2674 EPETRA_CHK_ERR(x_tmp->Import(*x_blocked,*Importer(), Insert)); // x_tmp will have all the values we need
2675 xp = (double*)x_tmp->Values();
2676 }
2677 }
2678
2679 for (int i=0; i < NumMyBlockRows_; i++) {
2680 int NumEntries = *NumBlockEntriesPerRow++;
2681 int * BlockRowIndices = *Indices++;
2682 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2683 int xoff = *RowFirstPointInElementList++;
2684 int RowDim = *RowElementSizeList++;
2685 if (DoRows) {
2686 for (int ii=0; ii < NumEntries; ii++) {
2687 double * xptr = xp+xoff;
2688 double * A = BlockRowValues[ii]->A();
2689 int LDA = BlockRowValues[ii]->LDA();
2690 int BlockIndex = BlockRowIndices[ii];
2691 int ColDim = ColElementSizeList[BlockIndex];
2692 for (int j=0; j<ColDim; j++) {
2693 double * curEntry = A + j*LDA;
2694 for (int k=0; k<RowDim; k++)
2695 *curEntry++ *= xptr[k];
2696 }
2697 }
2698 }
2699 else {
2700 for (int ii=0; ii < NumEntries; ii++) {
2701 double * A = BlockRowValues[ii]->A();
2702 int LDA = BlockRowValues[ii]->LDA();
2703 int BlockIndex = BlockRowIndices[ii];
2704 int off = ColFirstPointInElementList[BlockIndex];
2705 int ColDim = ColElementSizeList[BlockIndex];
2706 double * curx = xp+off;
2707 for (int j=0; j<ColDim; j++) {
2708 double * curEntry = A + j*LDA;
2709 for (int k=0; k<RowDim; k++)
2710 *curEntry++ *= curx[j];
2711 }
2712 }
2713 }
2714 }
2715
2716 if (x_tmp!=0) delete x_tmp;
2717 NormOne_ = -1.0; // Reset Norm so it will be recomputed.
2718 NormInf_ = -1.0; // Reset Norm so it will be recomputed.
2719 NormFrob_ = -1.0;
2721
2722 EPETRA_CHK_ERR(ierr);
2723 return(0);
2724}
2725//=============================================================================
2727
2728#if 0
2729 //
2730 // Commenting this section out disables caching, ie.
2731 // causes the norm to be computed each time NormInf is called.
2732 // See bug #1151 for a full discussion.
2733 //
2734 double MinNorm ;
2735 Comm().MinAll( &NormInf_, &MinNorm, 1 ) ;
2736
2737 if( MinNorm >= 0.0)
2738 // if( NormInf_ >= 0.0)
2739 return(NormInf_);
2740#endif
2741
2742 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
2743
2744 int MaxRowDim_ = MaxRowDim();
2745 double * tempv = new double[MaxRowDim_];
2746
2747 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2748 int * ElementSize = ElementSizeList_;
2750
2751 double Local_NormInf = 0.0;
2752 for (int i=0; i < NumMyBlockRows_; i++) {
2753 int NumEntries = *NumBlockEntriesPerRow++ ;
2754 int RowDim = *ElementSize++;
2755 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2756 BlockRowNormInf(RowDim, NumEntries,
2757 BlockRowValues, tempv);
2758 for (int j=0; j < RowDim; j++) Local_NormInf = EPETRA_MAX(Local_NormInf, tempv[j]);
2759 }
2760 Comm().MaxAll(&Local_NormInf, &NormInf_, 1);
2761 delete [] tempv;
2763 return(NormInf_);
2764}
2765//=============================================================================
2766void Epetra_VbrMatrix::BlockRowNormInf(int RowDim, int NumEntries,
2768 double * Y) const
2769{
2770 int i, j, k;
2771 for (k=0; k<RowDim; k++) Y[k] = 0.0;
2772
2773 for (i=0; i < NumEntries; i++) {
2774 double * A = As[i]->A();
2775 int LDA = As[i]->LDA();
2776 int ColDim = As[i]->N();
2777 for (j=0; j<ColDim; j++) {
2778 for (k=0; k<RowDim; k++) Y[k] += std::abs(A[k]);
2779 A += LDA;
2780 }
2781 }
2782 return;
2783}
2784//=============================================================================
2786
2787#if 0
2788 //
2789 // Commenting this section out disables caching, ie.
2790 // causes the norm to be computed each time NormOne is called.
2791 // See bug #1151 for a full discussion.
2792 //
2793 double MinNorm ;
2794 Comm().MinAll( &NormOne_, &MinNorm, 1 ) ;
2795
2796 if( MinNorm >= 0.0)
2797 return(NormOne_);
2798#endif
2799
2800 int * ColFirstPointInElementList = FirstPointInElementList_;
2801 if (Importer()!=0) {
2802 ColFirstPointInElementList = ColMap().FirstPointInElementList();
2803 }
2804
2805 Epetra_Vector * x = new Epetra_Vector(RowMap()); // Need temp vector for column sums
2806
2807 double * xp = (double*)x->Values();
2808 Epetra_MultiVector * x_tmp = 0;
2809
2810 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
2811 if (Importer()!=0) {
2812 x_tmp = new Epetra_Vector(ColMap()); // Create temporary import vector if needed
2813 xp = (double*)x_tmp->Values();
2814 }
2815
2816 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2817 int * ElementSize = ElementSizeList_;
2818 int ** Indices = Indices_;
2820
2821 for (int i=0; i < NumMyBlockRows_; i++) {
2822 int NumEntries = *NumBlockEntriesPerRow++;
2823 int RowDim = *ElementSize++;
2824 int * BlockRowIndices = *Indices++;
2825 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2826 BlockRowNormOne(RowDim, NumEntries, BlockRowIndices,
2827 BlockRowValues, ColFirstPointInElementList, xp);
2828 }
2829 if (Importer()!=0) {
2830 x->PutScalar(0.0);
2831 EPETRA_CHK_ERR(x->Export(*x_tmp, *Importer(), Add));
2832 } // Fill x with Values from temp vector
2833 x->MaxValue(&NormOne_); // Find max
2834 if (x_tmp!=0) delete x_tmp;
2835 delete x;
2837 return(NormOne_);
2838}
2839//=============================================================================
2841
2842#if 0
2843 //
2844 // Commenting this section out disables caching, ie.
2845 // causes the norm to be computed each time NormOne is called.
2846 // See bug #1151 for a full discussion.
2847 //
2848 double MinNorm ;
2849 Comm().MinAll( &NormFrob_, &MinNorm, 1 ) ;
2850
2851 if( MinNorm >= 0.0)
2852 return(NormFrob_);
2853#endif
2854
2855 int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
2856 int * ElementSize = ElementSizeList_;
2858
2859 double local_sum = 0.0;
2860
2861 for (int i=0; i < NumMyBlockRows_; i++) {
2862 int NumEntries = *NumBlockEntriesPerRow++;
2863 int RowDim = *ElementSize++;
2864 Epetra_SerialDenseMatrix** BlockRowValues = *Entries++;
2865
2866 for(int ii=0; ii<NumEntries; ++ii) {
2867 double* A = BlockRowValues[ii]->A();
2868 int LDA = BlockRowValues[ii]->LDA();
2869 int colDim = BlockRowValues[ii]->N();
2870 for(int j=0; j<colDim; ++j) {
2871 for(int k=0; k<RowDim; ++k) {
2872 local_sum += A[k]*A[k];
2873 }
2874 A += LDA;
2875 }
2876 }
2877
2878//The following commented-out code performs the calculation by running
2879//all the way across each point-entry row before going to the next.
2880//This is the order in which the CrsMatrix frobenius norm is calculated,
2881//but for a VbrMatrix it is faster to run the block-entries in
2882//column-major order (which is how they're stored), as the above code
2883//does.
2884//But if the same matrix is stored in both Vbr and Crs form, and you wish
2885//to compare the frobenius norms, it is necessary to run through the
2886//coefficients in the same order to avoid round-off differences.
2887//
2888// for(int k=0; k<RowDim; ++k) {
2889// for(int ii=0; ii<NumEntries; ++ii) {
2890// double* A = BlockRowValues[ii]->A();
2891// int LDA = BlockRowValues[ii]->LDA();
2892// int colDim = BlockRowValues[ii]->N();
2893// for(int j=0; j<colDim; ++j) {
2894// double Ak = A[k+j*LDA];
2895// local_sum += Ak*Ak;
2896// }
2897// }
2898// }
2899
2900 }
2901
2902 double global_sum = 0.0;
2903 Comm().SumAll(&local_sum, &global_sum, 1);
2904
2905 NormFrob_ = std::sqrt(global_sum);
2906
2908
2909 return(NormFrob_);
2910}
2911//=============================================================================
2912void Epetra_VbrMatrix::BlockRowNormOne(int RowDim, int NumEntries, int * BlockRowIndices,
2914 int * ColFirstPointInElementList, double * x) const {
2915 int i, j, k;
2916
2917 for (i=0; i < NumEntries; i++) {
2918 double * A = As[i]->A();
2919 int LDA = As[i]->LDA();
2920 int ColDim = As[i]->N();
2921 double * curx = x + ColFirstPointInElementList[BlockRowIndices[i]];
2922 for (j=0; j<ColDim; j++) {
2923 for (k=0; k<RowDim; k++) curx[j] += std::abs(A[k]);
2924 A += LDA;
2925 }
2926 }
2927 return;
2928}
2929//=========================================================================
2931 const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
2932 if (!A.Graph().GlobalConstantsComputed()) EPETRA_CHK_ERR(-1); // Must have global constants to proceed
2933 return(0);
2934}
2935//=========================================================================
2937 int NumSameIDs,
2938 int NumPermuteIDs,
2939 int * PermuteToLIDs,
2940 int *PermuteFromLIDs,
2941 const Epetra_OffsetIndex * Indexor,
2942 Epetra_CombineMode CombineMode)
2943{
2944 (void)Indexor;
2945 const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
2946 int i, j;
2947
2948 int BlockRow, NumBlockEntries;
2949 int * BlockIndices;
2950 int RowDim;
2951 Epetra_SerialDenseMatrix ** Entries;
2952 int FromBlockRow, ToBlockRow;
2953
2954 // Do copy first
2955 if (NumSameIDs>0) {
2956 int maxNumBlockEntries = A.MaxNumBlockEntries();
2957 BlockIndices = new int[maxNumBlockEntries]; // Need some temporary space
2958
2959
2960 for (i=0; i<NumSameIDs; i++) {
2961 BlockRow = GRID64(i);// FIXME long long
2963 maxNumBlockEntries,
2964 RowDim, NumBlockEntries,
2965 BlockIndices, Entries)); // Set pointers
2966 // Place into target matrix. Depends on Epetra_DataAccess copy/view and static/dynamic graph.
2967 if (StaticGraph() || IndicesAreLocal()) {
2968 if (CombineMode==Epetra_AddLocalAlso) {
2969 EPETRA_CHK_ERR(BeginSumIntoGlobalValues(BlockRow, NumBlockEntries,
2970 BlockIndices));
2971 }
2972 else {
2973 EPETRA_CHK_ERR(BeginReplaceGlobalValues(BlockRow, NumBlockEntries,
2974 BlockIndices));
2975 }
2976 }
2977 else {
2978 if (CombineMode==Epetra_AddLocalAlso) {
2979 EPETRA_CHK_ERR(BeginSumIntoGlobalValues(BlockRow, NumBlockEntries, BlockIndices));
2980 }
2981 else {
2982 EPETRA_CHK_ERR(BeginInsertGlobalValues(BlockRow, NumBlockEntries, BlockIndices));
2983 }
2984 }
2985 // Insert block entries one-at-a-time
2986 for (j=0; j<NumBlockEntries; j++) SubmitBlockEntry(Entries[j]->A(),
2987 Entries[j]->LDA(),
2988 RowDim, Entries[j]->N());
2989 EndSubmitEntries(); // Complete this block row
2990 }
2991 delete [] BlockIndices;
2992 }
2993
2994 // Do local permutation next
2995 if (NumPermuteIDs>0) {
2996 int maxNumBlockEntries = A.MaxNumBlockEntries();
2997 BlockIndices = new int[maxNumBlockEntries]; // Need some temporary space
2998
2999 for (i=0; i<NumPermuteIDs; i++) {
3000 FromBlockRow = A.GRID64(PermuteFromLIDs[i]);// FIXME long long
3001 ToBlockRow = GRID64(PermuteToLIDs[i]);// FIXME long long
3002 EPETRA_CHK_ERR(A.ExtractGlobalBlockRowPointers(FromBlockRow, maxNumBlockEntries, RowDim, NumBlockEntries,
3003 BlockIndices, Entries)); // Set pointers
3004 // Place into target matrix. Depends on Epetra_DataAccess copy/view and static/dynamic graph.
3005 if (StaticGraph() || IndicesAreLocal()) {
3006 if (CombineMode==Epetra_AddLocalAlso) {
3007 EPETRA_CHK_ERR(BeginSumIntoGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3008 }
3009 else {
3010 EPETRA_CHK_ERR(BeginReplaceGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3011 }
3012 }
3013 else {
3014 if (CombineMode==Epetra_AddLocalAlso) {
3015 EPETRA_CHK_ERR(BeginSumIntoGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3016 }
3017 else {
3018 EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3019 }
3020 }
3021 // Insert block entries one-at-a-time
3022 for (j=0; j<NumBlockEntries; j++) SubmitBlockEntry(Entries[j]->A(),
3023 Entries[j]->LDA(), RowDim, Entries[j]->N());
3024 EndSubmitEntries(); // Complete this block row
3025 }
3026 delete [] BlockIndices;
3027 }
3028
3029 return(0);
3030}
3031
3032//=========================================================================
3034 int NumExportIDs,
3035 int * ExportLIDs,
3036 int & LenExports,
3037 char * & Exports,
3038 int & SizeOfPacket,
3039 int * Sizes,
3040 bool & VarSizes,
3041 Epetra_Distributor & Distor)
3042{
3043 (void)LenExports;
3044 (void)Sizes;
3045 (void)VarSizes;
3046 (void)Distor;
3047 const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
3048
3049 double * DoubleExports = 0;
3050 int globalMaxNumNonzeros = A.GlobalMaxNumNonzeros();
3051 int globalMaxNumBlockEntries = A.GlobalMaxNumBlockEntries();
3052 // Will have globalMaxNumEntries doubles, globalMaxNumEntries +2 ints, pack them interleaved
3053 int DoublePacketSize = globalMaxNumNonzeros +
3054 (((2*globalMaxNumBlockEntries+3)+(int)sizeof(int)-1)*(int)sizeof(int))/(int)sizeof(double);
3055 SizeOfPacket = DoublePacketSize * (int)sizeof(double);
3056
3057 if (DoublePacketSize*NumExportIDs>LenExports_) {
3058 if (LenExports_>0) delete [] Exports_;
3059 LenExports_ = DoublePacketSize*NumExportIDs;
3060 DoubleExports = new double[LenExports_];
3061 Exports_ = (char *) DoubleExports;
3062 }
3063
3064 if (NumExportIDs<=0) return(0); // All done if nothing to pack
3065
3066 int i, j;
3067
3068 int NumBlockEntries;
3069 int * BlockIndices;
3070 int RowDim, * ColDims;
3071 double * Entries;
3072 int FromBlockRow;
3073 double * valptr, * dintptr;
3074 int * intptr;
3075
3076 // Each segment of IntExports will be filled by a packed row of information for each row as follows:
3077 // 1st int: GRID of row where GRID is the global row ID for the source matrix
3078 // next int: RowDim of Block Row
3079 // next int: NumBlockEntries, Number of indices in row.
3080 // next NumBlockEntries: The actual indices for the row.
3081 // Any remaining space (of length globalMaxNumNonzeros - NumBlockEntries ints) will be wasted but we need fixed
3082 // sized segments for current communication routines.
3083
3084 // Each segment of Exports will be filled with values.
3085
3086 valptr = (double *) Exports;
3087 dintptr = valptr + globalMaxNumNonzeros;
3088 intptr = (int *) dintptr;
3089 bool NoSumInto = false;
3090 for (i=0; i<NumExportIDs; i++) {
3091 FromBlockRow = A.GRID64(ExportLIDs[i]);// FIXME long long
3092 BlockIndices = intptr + 3;
3093 ColDims = BlockIndices + globalMaxNumBlockEntries;
3094 EPETRA_CHK_ERR(A.BeginExtractGlobalBlockRowCopy(FromBlockRow, globalMaxNumBlockEntries, RowDim,
3095 NumBlockEntries, BlockIndices, ColDims));
3096 // Now extract each block entry into send buffer
3097 Entries = valptr;
3098 for (j=0; j<NumBlockEntries; j++) {
3099 int SizeOfValues = RowDim*ColDims[j];
3100 A.ExtractEntryCopy(SizeOfValues, Entries, RowDim, NoSumInto);
3101 Entries += SizeOfValues;
3102 }
3103 // Fill first three slots of intptr with info
3104 intptr[0] = FromBlockRow;
3105 intptr[1] = RowDim;
3106 intptr[2] = NumBlockEntries;
3107 valptr += DoublePacketSize; // Point to next segment
3108 dintptr = valptr + globalMaxNumNonzeros;
3109 intptr = (int *) dintptr;
3110 }
3111
3112 return(0);
3113}
3114
3115//=========================================================================
3117 int NumImportIDs,
3118 int * ImportLIDs,
3119 int LenImports,
3120 char * Imports,
3121 int & SizeOfPacket,
3122 Epetra_Distributor & Distor,
3123 Epetra_CombineMode CombineMode,
3124 const Epetra_OffsetIndex * Indexor)
3125{
3126 (void)LenImports;
3127 (void)SizeOfPacket;
3128 (void)Distor;
3129 (void)Indexor;
3130 if (NumImportIDs<=0) return(0);
3131
3132 if( CombineMode != Add
3133 && CombineMode != Zero
3134 && CombineMode != Insert )
3135 EPETRA_CHK_ERR(-1); // CombineMode not supported, default to mode Zero
3136
3137 const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
3138 int NumBlockEntries;
3139 int * BlockIndices;
3140 int RowDim, * ColDims;
3141 double * values;
3142 int ToBlockRow;
3143 int i, j;
3144
3145 double * valptr, *dintptr;
3146 int * intptr;
3147 int globalMaxNumNonzeros = A.GlobalMaxNumNonzeros();
3148 int globalMaxNumBlockEntries = A.GlobalMaxNumBlockEntries();
3149 // Will have globalMaxNumEntries doubles, globalMaxNumEntries +2 ints, pack them interleaved
3150 int DoublePacketSize = globalMaxNumNonzeros +
3151 (((2*globalMaxNumBlockEntries+3)+(int)sizeof(int)-1)*(int)sizeof(int))/(int)sizeof(double);
3152 // Unpack it...
3153
3154
3155 // Each segment of IntImports will be filled by a packed row of information for each row as follows:
3156 // 1st int: GRID of row where GRID is the global row ID for the source matrix
3157 // next int: NumBlockEntries, Number of indices in row.
3158 // next int: RowDim of Block Row
3159 // next NumBlockEntries: The actual indices for the row.
3160 // Any remaining space (of length globalMaxNumNonzeros - NumBlockEntries ints) will be
3161 // wasted but we need fixed sized segments for current communication routines.
3162
3163 valptr = (double *) Imports;
3164 dintptr = valptr + globalMaxNumNonzeros;
3165 intptr = (int *) dintptr;
3166
3167 for (i=0; i<NumImportIDs; i++) {
3168 ToBlockRow = GRID64(ImportLIDs[i]);// FIXME long long
3169 assert((intptr[0])==ToBlockRow); // Sanity check
3170 RowDim = RowMap().ElementSize(ImportLIDs[i]);
3171 assert((intptr[1])==RowDim); // Sanity check
3172 NumBlockEntries = intptr[2];
3173 BlockIndices = intptr + 3;
3174 ColDims = BlockIndices + globalMaxNumBlockEntries;
3175 if (CombineMode==Add) {
3176 if (StaticGraph() || IndicesAreLocal()) {
3177 // Replace any current values
3178 EPETRA_CHK_ERR(BeginSumIntoGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3179 }
3180 else {
3181 // Insert values
3182 EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3183 }
3184 }
3185 else if (CombineMode==Insert) {
3186 if (StaticGraph() || IndicesAreLocal()) {
3187 // Replace any current values
3188 EPETRA_CHK_ERR(BeginReplaceGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3189 }
3190 else {
3191 // Insert values
3192 EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
3193 }
3194 }
3195 // Now extract each block entry into send buffer
3196 values = valptr;
3197 for (j=0; j<NumBlockEntries; j++) {
3198 int LDA = RowDim;
3199 int ColDim = ColDims[j];
3200 SubmitBlockEntry(values, LDA, RowDim, ColDim);
3201 values += (LDA*ColDim);
3202 }
3203 EndSubmitEntries(); // Done with this block row
3204 valptr += DoublePacketSize; // Point to next segment
3205 dintptr = valptr + globalMaxNumNonzeros;
3206 intptr = (int *) dintptr;
3207 }
3208
3209 return(0);
3210}
3211//=========================================================================
3213
3214 if (HavePointObjects_) return(0); // Already done
3215
3216 // Generate a point-wise compatible row map
3218
3219 // For each of the next maps, first check if the corresponding block map
3220 // is the same as the block row map for this matrix. If so, then we simply
3221 // copy the pointer. Otherwise we create a new point map.
3222
3223 // This check can save storage and time since it will often be the case that the
3224 // domain, range and row block maps will be the same. Also, in the serial case,
3225 // the column block map will also often be the same as the row block map.
3226
3227 if (ColMap().SameAs(RowMap())) {
3229 }
3230 else {
3232 }
3233
3234 if (DomainMap().SameAs(RowMap())) {
3236 }
3237 else {
3239 }
3240 if (RangeMap().SameAs(RowMap())) {
3242 }
3243 else {
3245 }
3246
3247 // Finally generate Importer that will migrate needed domain elements to the column map
3248 // layout.
3250
3251 HavePointObjects_ = true;
3252 return(0);
3253}
3254//=========================================================================
3256 Epetra_Map * & PointMap) const
3257{
3258 // Generate an Epetra_Map that has the same number and distribution of points
3259 // as the input Epetra_BlockMap object. The global IDs for the output PointMap
3260 // are computed by using the MaxElementSize of the BlockMap. For variable block
3261 // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.
3262
3263 int MaxElementSize = BlockMap.MaxElementSize();
3264 int PtNumMyElements = BlockMap.NumMyPoints();
3265 int * PtMyGlobalElements = 0;
3266 if (PtNumMyElements>0) PtMyGlobalElements = new int[PtNumMyElements];
3267
3268 int NumMyElements = BlockMap.NumMyElements();
3269
3270 int curID = 0;
3271 for (int i=0; i<NumMyElements; i++) {
3272 int StartID = BlockMap.GID64(i)*MaxElementSize; // CJ TODO FIXME long long
3273 int ElementSize = BlockMap.ElementSize(i);
3274 for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
3275 }
3276 assert(curID==PtNumMyElements); // Sanity test
3277
3278 PointMap = new Epetra_Map(-1, PtNumMyElements, PtMyGlobalElements, BlockMap.IndexBase(), BlockMap.Comm());// FIXME long long
3279
3280 if (PtNumMyElements>0) delete [] PtMyGlobalElements;
3281
3282 if (!BlockMap.PointSameAs(*PointMap)) {EPETRA_CHK_ERR(-1);} // Maps not compatible
3283 return(0);
3284}
3285
3286//=========================================================================
3288 const Epetra_MultiVector& Y) const
3289{
3290 if (OperatorX_!=0)
3291 if (OperatorX_->NumVectors()!=X.NumVectors()) {delete OperatorX_; OperatorX_ = 0; delete OperatorY_; OperatorY_=0;}
3292 if (OperatorX_==0) {
3293 if (!X.Map().PointSameAs(DomainMap())) EPETRA_CHK_ERR(-1); // X not point-wise compatible with the block domain map
3294 if (!Y.Map().PointSameAs(RangeMap())) EPETRA_CHK_ERR(-2); // Y not point-wise compatible with the block col map
3297 }
3298 else {
3301 }
3302 return(0);
3303}
3304//=========================================================================
3306 Epetra_MultiVector& Y) const
3307{
3309 EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3311 }
3312 else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
3313 EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3315 }
3316 return(0);
3317}
3318//=========================================================================
3320 Epetra_MultiVector& Y) const
3321{
3323 EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3325 }
3326 else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
3327 EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3329 }
3330 return(0);
3331}
3332//=========================================================================
3334 Epetra_MultiVector& Y) const
3335{
3336 if (!TransA) {
3337 EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3339 }
3340 else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
3341 EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3343 }
3344 return(0);
3345}
3346//=========================================================================
3347int Epetra_VbrMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X,
3348 Epetra_MultiVector& Y) const
3349{
3350 if (!Trans) {
3351 EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3352 EPETRA_CHK_ERR(DoSolve(Upper, Trans, UnitDiagonal, *OperatorX_, *OperatorY_));
3353 }
3354 else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
3355 EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
3356 EPETRA_CHK_ERR(DoSolve(Upper, Trans, UnitDiagonal, *OperatorY_, *OperatorX_));
3357 }
3358 return(0);
3359}
3360//=========================================================================
3361void Epetra_VbrMatrix::Print(std::ostream& os) const {
3362 int MyPID = RowMap().Comm().MyPID();
3363 int NumProc = RowMap().Comm().NumProc();
3364
3365 for (int iproc=0; iproc < NumProc; iproc++) {
3366 if (MyPID==iproc) {
3367 if (MyPID==0) {
3368 os << "\nNumber of Global Block Rows = "; os << NumGlobalBlockRows64(); os << std::endl;
3369 os << "Number of Global Block Cols = "; os << NumGlobalBlockCols64(); os << std::endl;
3370 os << "Number of Global Block Diags = "; os << NumGlobalBlockDiagonals64(); os << std::endl;
3371 os << "Number of Global Blk Entries = "; os << NumGlobalBlockEntries64(); os << std::endl;
3372 os << "Global Max Num Block Entries = "; os << GlobalMaxNumBlockEntries(); os << std::endl;
3373 os << "\nNumber of Global Rows = "; os << NumGlobalRows64(); os << std::endl;
3374 os << "Number of Global Cols = "; os << NumGlobalCols64(); os << std::endl;
3375 os << "Number of Global Diagonals = "; os << NumGlobalDiagonals64(); os << std::endl;
3376 os << "Number of Global Nonzeros = "; os << NumGlobalNonzeros64(); os << std::endl;
3377 os << "Global Maximum Num Entries = "; os << GlobalMaxNumNonzeros(); os << std::endl;
3378 if (LowerTriangular()) { os << " ** Matrix is Lower Triangular **"; os << std::endl; }
3379 if (UpperTriangular()) { os << " ** Matrix is Upper Triangular **"; os << std::endl; }
3380 if (NoDiagonal()) { os << " ** Matrix has no diagonal **"; os << std::endl; os << std::endl; }
3381 }
3382
3383 os << "\nNumber of My Block Rows = "; os << NumMyBlockRows(); os << std::endl;
3384 os << "Number of My Block Cols = "; os << NumMyBlockCols(); os << std::endl;
3385 os << "Number of My Block Diags = "; os << NumMyBlockDiagonals(); os << std::endl;
3386 os << "Number of My Blk Entries = "; os << NumMyBlockEntries(); os << std::endl;
3387 os << "My Max Num Block Entries = "; os << MaxNumBlockEntries(); os << std::endl;
3388 os << "\nNumber of My Rows = "; os << NumMyRows(); os << std::endl;
3389 os << "Number of My Cols = "; os << NumMyCols(); os << std::endl;
3390 os << "Number of My Diagonals = "; os << NumMyDiagonals(); os << std::endl;
3391 os << "Number of My Nonzeros = "; os << NumMyNonzeros(); os << std::endl;
3392 os << "My Maximum Num Entries = "; os << MaxNumBlockEntries(); os << std::endl; os << std::endl;
3393
3394 os << std::flush;
3395
3396 }
3397 // Do a few global ops to give I/O a chance to complete
3398 Comm().Barrier();
3399 Comm().Barrier();
3400 Comm().Barrier();
3401 }
3402
3403 {for (int iproc=0; iproc < NumProc; iproc++) {
3404 if (MyPID==iproc) {
3405 int NumBlockRows1 = NumMyBlockRows();
3406 int MaxNumBlockEntries1 = MaxNumBlockEntries();
3407 int * BlockIndices1 = new int[MaxNumBlockEntries1];
3408 Epetra_SerialDenseMatrix ** Entries1;
3409 int RowDim1, NumBlockEntries1;
3410 int i, j;
3411
3412 if (MyPID==0) {
3413 os.width(8);
3414 os << " Processor ";
3415 os.width(10);
3416 os << " Block Row Index ";
3417 os.width(10);
3418 os << " Block Col Index \n";
3419 os.width(20);
3420 os << " values ";
3421 os << std::endl;
3422 }
3423 for (i=0; i<NumBlockRows1; i++) {
3424 int BlockRow1 = GRID64(i); // Get global row number// FIXME long long
3425 ExtractGlobalBlockRowPointers(BlockRow1, MaxNumBlockEntries1, RowDim1,
3426 NumBlockEntries1, BlockIndices1,
3427 Entries1);
3428
3429 for (j = 0; j < NumBlockEntries1 ; j++) {
3430 os.width(8);
3431 os << MyPID ; os << " ";
3432 os.width(10);
3433 os << BlockRow1 ; os << " ";
3434 os.width(10);
3435 os << BlockIndices1[j]; os << " " << std::endl;
3436 os.width(20);
3437
3438 if (Entries1[j] == 0) {
3439 os << "Block Entry == NULL"<< std::endl;
3440 continue;
3441 }
3442
3443 Epetra_SerialDenseMatrix entry_view(View, Entries1[j]->A(), Entries1[j]->LDA(),
3444 RowDim1, Entries1[j]->N());
3445 os << entry_view; os << " ";
3446 os << std::endl;
3447 }
3448 }
3449
3450 delete [] BlockIndices1;
3451
3452 os << std::flush;
3453
3454 }
3455 // Do a few global ops to give I/O a chance to complete
3456 Comm().Barrier();
3457 Comm().Barrier();
3458 Comm().Barrier();
3459 }}
3460
3461 return;
3462}
3463
3464#endif // EPETRA_NO_32BIT_GLOBAL_INDICES
Epetra_CombineMode
@ Insert
@ Zero
@ Epetra_AddLocalAlso
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
const double Epetra_MinDouble
const double Epetra_MaxDouble
Epetra_DataAccess
@ View
@ Copy
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
void GEMV(const char TRANS, const int M, const int N, const float ALPHA, const float *A, const int LDA, const float *X, const float BETA, float *Y, const int INCX=1, const int INCY=1) const
Epetra_BLAS matrix-vector multiply function (SGEMV)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MaxElementSize() const
Maximum element size across all processors.
bool PointSameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map have identical point-wise structure.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor.
int IndexBase() const
Index base for this map.
long long GID64(int LID) const
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
Returns the LID of the element that contains the given local PointID, and the Offset of the point in ...
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
bool MyLID(int lid) const
Returns true if the LID passed in belongs to the calling processor in this map, otherwise returns fal...
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
virtual int MyPID() const =0
Return my process ID.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
Epetra_CompObject: Functionality and data that is common to all computational classes.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false.
int * NumIndicesPerRow() const
bool FindGlobalIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
int MakeIndicesLocal(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
void SetIndicesAreLocal(bool Flag)
int ** Indices() const
bool HaveColMap() const
Returns true if we have a well-defined ColMap, and returns false otherwise.
bool GlobalConstantsComputed() const
int ExtractGlobalRowView(int GlobalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified global row of the graph.
void SetSorted(bool Flag)
int ExtractMyRowCopy(int LocalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified local row of the graph. Put into storage allocated by calli...
int RemoveRedundantIndices()
Removes any redundant column indices in the rows of the graph.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int InsertIndices(int Row, int NumIndices, int *Indices)
int NumMyEntries() const
Returns the number of entries on this processor.
bool FindMyIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
void SetIndicesAreGlobal(bool Flag)
int ExtractGlobalRowCopy(int GlobalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified global row of the graph. Put into storage allocated by call...
int * NumAllocatedIndicesPerRow() const
const Epetra_BlockMap & ColMap() const
Returns the Column Map associated with this graph.
int NumMyNonzeros() const
Returns the number of indices in the local graph.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
Definition: Epetra_Import.h:63
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumVectors() const
Returns the number of vectors in the multi-vector.
int ResetView(double **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
double * Values() const
Get pointer to MultiVector values.
double ** Pointers() const
Get pointer to individual vector pointers.
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
virtual int ColDim() const
Returns the column dimension of operator.
void CopyMat(const double *Source, int Source_LDA, int NumRows, int NumCols, double *Target, int Target_LDA, bool add=false)
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
virtual int RowDim() const
Returns the row dimension of operator.
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
const Epetra_Export * Exporter() const
Returns the Epetra_Export object that contains the export operations for distributed operations.
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false.
void BlockRowNormInf(int RowDim, int NumEntries, Epetra_SerialDenseMatrix **As, double *Y) const
int UpdateOperatorXY(const Epetra_MultiVector &X, const Epetra_MultiVector &Y) const
int NumMyBlockEntries() const
Returns the number of nonzero block entries in the calling processor's portion of the matrix.
int ExtractEntryView(Epetra_SerialDenseMatrix *&entry) const
Returns a pointer to the current block entry.
int InvRowSums(Epetra_Vector &x) const
Computes the sum of absolute values of the rows of the Epetra_VbrMatrix, results returned in x.
int BeginExtractGlobalBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices) const
Initiates a view of the specified global row, only works if matrix indices are in global mode.
virtual void Print(std::ostream &os) const
Print method.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
int DoMultiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
int ExtractMyBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Copy the block indices into user-provided array, set pointers for rest of data for specified local bl...
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
long long NumGlobalCols64() const
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
int ExtractBlockDiagonalEntryCopy(int SizeOfValues, double *Values, int LDA, bool SumInto) const
Extract a copy of a block diagonal entry from the matrix.
int PutScalar(double ScalarConstant)
Initialize all values in graph of the matrix with constant value.
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine.
Epetra_VbrMatrix & operator=(const Epetra_VbrMatrix &src)
long long NumGlobalBlockRows64() const
Epetra_Import * RowMatrixImporter_
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
int Multiply1(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_Vector x in y.
Epetra_MultiVector * OperatorX_
int BeginExtractMyBlockRowCopy(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, int *ColDims) const
Initiates a copy of the specified local row in user-provided arrays.
long long GRID64(int LRID_in) const
int BeginInsertMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given local row of the matrix, values are inserted via ...
int NumMyNonzeros() const
Returns the number of nonzero entriesowned by the calling processor .
void FastBlockRowMultiply(bool TransA, int RowDim, int NumEntries, int *BlockIndices, int RowOff, int *FirstPointInElementList, int *ElementSizeList, Epetra_SerialDenseMatrix **As, double **X, double **Y, int NumVectors) const
bool NoDiagonal() const
If matrix has no diagonal entries based on global row/column index comparisons, this query returns tr...
int BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate summing into current values with this list of entries for a given global row of the matrix,...
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not.
int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
int ExtractGlobalBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Copy the block indices into user-provided array, set pointers for rest of data for specified global b...
long long NumGlobalBlockDiagonals64() const
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
int BeginExtractBlockDiagonalView(int &NumBlockDiagonalEntries, int *&RowColDims) const
Initiates a view of the block diagonal entries.
int MaxRowDim() const
Returns the maximum row dimension of all block entries on this processor.
int ReplaceMatDiag(double *A, int LDA, int NumRows, int NumCols, double *Diagonal)
int SortEntries()
Sort column entries, row-by-row, in ascending order.
int GlobalMaxNumBlockEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
int ExtractEntryCopy(int SizeOfValues, double *Values, int LDA, bool SumInto) const
Extract a copy of an entry from the block row specified by one of the BeginExtract routines.
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the with those in the user-provided vector.
Epetra_MultiVector * ImportVector_
int BeginExtractMyBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices) const
Initiates a view of the specified local row, only works if matrix indices are in local mode.
virtual ~Epetra_VbrMatrix()
Epetra_VbrMatrix Destructor.
int MergeRedundantEntries()
Add entries that have the same column index. Remove redundant entries from list.
int BlockMap2PointMap(const Epetra_BlockMap &BlockMap, Epetra_Map *&PointMap) const
bool NoRedundancies() const
If MergeRedundantEntries() has been called, this query returns true, otherwise it returns false.
bool UseTranspose() const
Returns the current UseTranspose setting.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int CopyMatDiag(double *A, int LDA, int NumRows, int NumCols, double *Diagonal) const
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given global row of the matrix, values are inserted via...
int BeginReplaceGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate replacement of current values with this list of entries for a given global row of the matrix...
const Epetra_CrsGraph & Graph() const
Returns a pointer to the Epetra_CrsGraph object associated with this matrix.
Epetra_CrsGraph * Graph_
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
const Epetra_Comm & Comm() const
Fills a matrix with rows from a source matrix based on the specified importer.
int NumMyBlockDiagonals() const
Returns the number of local nonzero block diagonal entries, based on global row/column index comparis...
int BeginExtractBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, bool IndicesAreLocal) const
void BlockRowMultiply(bool TransA, int RowDim, int NumEntries, int *BlockIndices, int RowOff, int *FirstPointInElementList, int *ElementSizeList, double Alpha, Epetra_SerialDenseMatrix **As, double **X, double Beta, double **Y, int NumVectors) const
long long NumGlobalBlockEntries64() const
int MaxNumBlockEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
int SetAllocated(bool Flag)
bool Sorted() const
If SortEntries() has been called, this query returns true, otherwise it returns false.
void BlockRowNormOne(int RowDim, int NumEntries, int *BlockRowIndices, Epetra_SerialDenseMatrix **As, int *ColFirstPointInElementList, double *x) const
const Epetra_Import * Importer() const
Returns the Epetra_Import object that contains the import operations for distributed operations.
int NumMyBlockCols() const
Returns the number of Block matrix columns owned by the calling processor.
int FillComplete()
Signal that data entry is complete, perform transformations to local index space.
const Epetra_BlockMap & ColMap() const
Returns the ColMap as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing Epetra_R...
Epetra_Map * RowMatrixRowMap_
int NumMyBlockRows() const
Returns the number of Block matrix rows owned by the calling processor.
int ExtractMyBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Initiates a view of the specified local row, only works if matrix indices are in local mode.
long long NumGlobalBlockCols64() const
int BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate replacement of current values with this list of entries for a given local row of the matrix,...
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_MultiVector X in Y.
double NormInf() const
Returns the infinity norm of the global matrix.
Epetra_MultiVector * ExportVector_
int NumMyCols() const
Returns the number of matrix columns owned by the calling processor.
int ExtractBlockDiagonalEntryView(double *&Values, int &LDA) const
Extract a view of a block diagonal entry from the matrix.
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true,...
int ExtractGlobalBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, Epetra_SerialDenseMatrix **&Values) const
Initiates a view of the specified global row, only works if matrix indices are in global mode.
int BeginInsertValues(int BlockRow, int NumBlockEntries, int *BlockIndices, bool IndicesAreLocal)
int BeginExtractBlockDiagonalCopy(int MaxNumBlockDiagonalEntries, int &NumBlockDiagonalEntries, int *RowColDims) const
Initiates a copy of the block diagonal entries to user-provided arrays.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
int DoSolve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Epetra_DataAccess CV_
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
int BeginExtractBlockRowCopy(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, int *ColDims, bool IndicesAreLocal) const
long long NumGlobalNonzeros64() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int ExtractBlockDimsCopy(int NumBlockEntries, int *ColDims) const
int BeginReplaceValues(int BlockRow, int NumBlockEntries, int *BlockIndices, bool IndicesAreLocal)
int BeginSumIntoMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate summing into current values with this list of entries for a given local row of the matrix,...
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the left with a Epetra_Vector x.
int ExtractBlockRowPointers(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, Epetra_SerialDenseMatrix **&Values, bool IndicesAreLocal) const
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
int DirectSubmitBlockEntry(int GlobalBlockRow, int GlobalBlockCol, const double *values, int LDA, int NumRows, int NumCols, bool sum_into)
Submit a block-entry directly into the matrix (without using a begin/end sequence)
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a solve using the Epetra_VbrMatrix on a Epetra_Vector x in y.
int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
Perform any packing or preparation required for call to DoTransfer().
int CopyMat(double *A, int LDA, int NumRows, int NumCols, double *B, int LDB, bool SumInto) const
const Epetra_Map & RowMatrixColMap() const
Returns the Epetra_Map object associated with columns of this matrix.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
Epetra_SerialDenseMatrix *** Entries_
int RightScale(const Epetra_Vector &x)
Scales the Epetra_VbrMatrix on the right with a Epetra_Vector x.
int OptimizeStorage()
Eliminates memory that is used for construction. Make consecutive row index sections contiguous.
double NormOne() const
Returns the one norm of the global matrix.
Epetra_Map * OperatorDomainMap_
Epetra_SerialDenseMatrix ** TempEntries_
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A).
int BeginSumIntoValues(int BlockRow, int NumBlockEntries, int *BlockIndices, bool IndicesAreLocal)
const Epetra_BlockMap & RangeMap() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Epetra_Map * OperatorRangeMap_
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
int GeneratePointObjects() const
long long NumGlobalDiagonals64() const
int BeginExtractGlobalBlockRowCopy(int BlockRow, int MaxNumBlockEntries, int &RowDim, int &NumBlockEntries, int *BlockIndices, int *ColDims) const
Initiates a copy of the specified global row in user-provided arrays.
Epetra_Map * RowMatrixColMap_
bool StaticGraph() const
long long GCID64(int LCID_in) const
int SetupForExtracts(int BlockRow, int &RowDim, int NumBlockEntries, bool ExtractView, bool IndicesAreLocal) const
const Epetra_Map & RowMatrixRowMap() const
Returns the EpetraMap object associated with the rows of this matrix.
long long NumGlobalRows64() const
Epetra_CombineMode CurSubmitMode_
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
int NumMyDiagonals() const
Returns the number of local nonzero diagonal entries, based on global row/column index comparisons.
int TransformToLocal()
Use FillComplete() instead.
Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap &RowMap, int *NumBlockEntriesPerRow)
Epetra_VbrMatrix constuctor with variable number of indices per row.
int * NumAllocatedBlockEntriesPerRow_
int InvColSums(Epetra_Vector &x) const
Computes the sum of absolute values of the columns of the Epetra_VbrMatrix, results returned in x.
int GlobalMaxNumNonzeros() const
Returns the maximum number of nonzero entries across all block rows on all processors.
const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
int InverseSums(bool DoRows, Epetra_Vector &x) const
int SetupForSubmits(int BlockRow, int NumBlockEntries, int *BlockIndices, bool IndicesAreLocal, Epetra_CombineMode SubmitMode)
Epetra_MultiVector * OperatorY_
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.