Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_FastCrsMatrix.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
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"
52#include <limits.h>
53
54//==============================================================================
55Epetra_FastCrsMatrix::Epetra_FastCrsMatrix(const Epetra_CrsMatrix & Matrix, bool UseFloats)
56 : CrsMatrix_(Matrix),
57 Values_(0),
58 NumMyRows_(Matrix.NumMyRows()),
59 NumMyNonzeros(Matrix.NumMyNonzeros()),
60 ImportVector_(0),
61 ExportVector_(0),
62 CV_(Copy)
63{
64 if (!CrsMatrix_.Filled()) throw CrsMatrix_.ReportError("Input matrix must have called FillComplete()", -1);
65 Allocate(UseFloats);
66}
67
68//==============================================================================
69Epetra_FastCrsMatrix::~Epetra_FastCrsMatrix(){
70
71 if (Values_!=0 && ValuesAllocated_) delete [] Values_;
72 if (FloatValues_!=0) delete [] FloatValues_;
73 if (Indices_!=0) delete [] Indices_;
74 if (ShortIndices_!=0) delete [] ShortIndices_;
75
76 if (ImportVector_!=0) delete ImportVector_;
77 ImportVector_=0;
78 if (ExportVector_!=0) delete ExportVector_;
79 ExportVector_=0;
80}
81
82//==============================================================================
83int Epetra_FastCrsMatrix::UpdateValues(const Epetra_CrsMatrix & Matrix) {
84}
85
86//==============================================================================
87int Epetra_FastCrsMatrix::Allocate(bool UseFloats) {
88
89 int i, j;
90
91 // First determine how to handle values. Three possibilities:
92 // 1) Values are contiguously stored in user matrix, UseFloats is false so we point to the values in
93 // the user matrix.
94 // 2) Values are not contiguously stored, UseFloats is false, so we copy values into a contiguous double array.
95 // 3) UseFloats is true so we create a single precision copy of matrix values.
96
97 if (CrsMatrix.IndicesAreContiguous() && !UseFloats) {
98 ValuesAllocated_ = false;
99 Values_ = CrsMatrix_[0]; // First value in the user's matrix
100 }
101 else if (!UseFloats) {
102 // Allocate Values array
103 Values_ = new double[NumMyNonzeros_];
104 double * ptr = Values_;
105 for (i=0; i<NumMyRows_; i++) {
106 int NumEntries = CrsMatrix_.NumMyEntries(i);
107 double * rowi = CrsMatrix_[i];
108 for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
109 }
110 }
111 else { // UseFloats == true
112 FloatValues_ = new float[NumMyNonzeros_];
113 float * ptr = FloatValues_;
114 for (i=0; i<NumMyRows_; i++) {
115 int NumEntries = CrsMatrix_.NumMyEntries(i);
116 double * rowi = CrsMatrix_[i];
117 for (j=0; j< NumEntries; j++) *ptr++ = (float) *rowi++; // convert and copy values
118 }
119 }
120
121 // Next determine how to handle integers. Two possibilities:
122 // 1) Local column range is within the range of unsigned short ints, so we copy the indices to an array of unsigned shorts.
123 // 2) Local column range is outside range of unsigned shorts and we copy to array of ints.
124 // In both cases we embed the nonzero count per row into the index array.
125
126 if (CrsMatrix_.NumMyCols()>USHRT_MAX) {
127 // Allocate Values array
128 Indices_ = new int[NumMyNonzeros_+NumMyRows_];
129 int * ptr = Indices_;
130 for (i=0; i<NumMyRows_; i++) {
131 int NumEntries = CrsMatrix_.NumMyEntries(i);
132 int * rowi = CrsMatrix_.Graph()[i];
133 *ptr++ = NumEntries; // Pack this value first
134 for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
135 }
136 }
137 else { // CrsMatrix_.NumMyCols()<=USHRT_MAX
138 ShortIndices_ = new unsigned short[NumMyNonzeros_+NumMyRows_];
139 unsigned short * ptr = ShortIndices_;
140 for (i=0; i<NumMyRows_; i++) {
141 int NumEntries = CrsMatrix_.NumMyEntries(i);
142 unsigned short * rowi = CrsMatrix_Graph()[i];
143 *ptr++ = NumEntries; // Pack this value first
144 for (j=0; j< NumEntries; j++) *ptr++ = (unsigned short) *rowi++; // convert and copy indices
145 }
146 }
147
148 return(0);
149}
150//=============================================================================
151int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
152 //
153 // This function forms the product y = A * x or y = A' * x
154 //
155
156 int i, j;
157 double * xp = (double*)x.Values();
158 double *yp = (double*)y.Values();
159 int NumMyCols_ = NumMyCols();
160
161
162 if (!TransA) {
163
164 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
165 if (Importer()!=0) {
166 if (ImportVector_!=0) {
167 if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
168 }
169 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
170 ImportVector_->Import(x, *Importer(), Insert);
171 xp = (double*)ImportVector_->Values();
172 }
173
174 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
175 if (Exporter()!=0) {
176 if (ExportVector_!=0) {
177 if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
178 }
179 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
180 yp = (double*)ExportVector_->Values();
181 }
182
183 // Do actual computation
184
185 for (i=0; i < NumMyRows_; i++) {
186 int NumEntries = *NumEntriesPerRow++;
187 int * RowIndices = *Indices++;
188 double * RowValues = *Values++;
189 double sum = 0.0;
190 for (j=0; j < NumEntries; j++) sum += RowValues[j] * xp[RowIndices[j]];
191
192 yp[i] = sum;
193
194 }
195 if (Exporter()!=0) y.Export(*ExportVector_, *Exporter(), Add); // Fill y with Values from export vector
196 }
197
198 else { // Transpose operation
199
200
201 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
202
203 if (Exporter()!=0) {
204 if (ExportVector_!=0) {
205 if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
206 }
207 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
208 ExportVector_->Import(x, *Exporter(), Insert);
209 xp = (double*)ExportVector_->Values();
210 }
211
212 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
213 if (Importer()!=0) {
214 if (ImportVector_!=0) {
215 if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
216 }
217 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
218 yp = (double*)ImportVector_->Values();
219 }
220
221 // Do actual computation
222
223 for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
224
225 for (i=0; i < NumMyRows_; i++) {
226 int NumEntries = *NumEntriesPerRow++;
227 int * RowIndices = *Indices++;
228 double * RowValues = *Values++;
229 for (j=0; j < NumEntries; j++) yp[RowIndices[j]] += RowValues[j] * xp[i];
230 }
231 if (Importer()!=0) y.Export(*ImportVector_, *Importer(), Add); // Fill y with Values from export vector
232 }
233
234 UpdateFlops(2*NumGlobalNonzeros64());
235 return(0);
236}
237
238//=============================================================================
239int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
240 //
241 // This function forms the product Y = A * Y or Y = A' * X
242 //
243 if (X.NumVectors()==1 && Y.NumVectors()==1) {
244 double * xp = (double *) X[0];
245 double * yp = (double *) Y[0];
246 Epetra_Vector x(View, X.Map(), xp);
247 Epetra_Vector y(View, Y.Map(), yp);
248 return(Multiply(TransA, x, y));
249 }
250 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
251
252 int i, j, k;
253 int * NumEntriesPerRow = NumEntriesPerRow_;
254 int ** Indices = Indices_;
255 double ** Values = Values_;
256
257 double **Xp = (double**)X.Pointers();
258 double **Yp = (double**)Y.Pointers();
259
260 int NumVectors = X.NumVectors();
261 int NumMyCols_ = NumMyCols();
262
263
264 // Need to better manage the Import and Export vectors:
265 // - Need accessor functions
266 // - Need to make the NumVector match (use a View to do this)
267 // - Need to look at RightScale and ColSum routines too.
268
269 if (!TransA) {
270
271 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
272 if (Importer()!=0) {
273 if (ImportVector_!=0) {
274 if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
275 }
276 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
277 ImportVector_->Import(X, *Importer(), Insert);
278 Xp = (double**)ImportVector_->Pointers();
279 }
280
281 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
282 if (Exporter()!=0) {
283 if (ExportVector_!=0) {
284 if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
285 }
286 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
287 Yp = (double**)ExportVector_->Pointers();
288 }
289
290 // Do actual computation
291
292 for (i=0; i < NumMyRows_; i++) {
293 int NumEntries = *NumEntriesPerRow++;
294 int * RowIndices = *Indices++;
295 double * RowValues = *Values++;
296 for (k=0; k<NumVectors; k++) {
297 double sum = 0.0;
298 for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
299 Yp[k][i] = sum;
300 }
301 }
302 if (Exporter()!=0) Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
303 }
304 else { // Transpose operation
305
306
307 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
308
309 if (Exporter()!=0) {
310 if (ExportVector_!=0) {
311 if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
312 }
313 if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
314 ExportVector_->Import(X, *Exporter(), Insert);
315 Xp = (double**)ExportVector_->Pointers();
316 }
317
318 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
319 if (Importer()!=0) {
320 if (ImportVector_!=0) {
321 if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
322 }
323 if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
324 Yp = (double**)ImportVector_->Pointers();
325 }
326
327 // Do actual computation
328
329
330
331 for (k=0; k<NumVectors; k++)
332 for (i=0; i < NumMyCols_; i++) Yp[k][i] = 0.0; // Initialize y for transpose multiply
333
334 for (i=0; i < NumMyRows_; i++) {
335 int NumEntries = *NumEntriesPerRow++;
336 int * RowIndices = *Indices++;
337 double * RowValues = *Values++;
338 for (k=0; k<NumVectors; k++) {
339 for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
340 }
341 }
342 if (Importer()!=0) Y.Export(*ImportVector_, *Importer(), Add); // Fill Y with Values from export vector
343 }
344
345 UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
346 return(0);
347}
348
349//=============================================================================
350int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
351 //
352 // This function find y such that Ly = x or Uy = x or the transpose cases.
353 //
354
355 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
356
357 if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
358 if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
359 if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
360 if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
361
362
363 int i, j, j0;
364 int * NumEntriesPerRow = NumEntriesPerRow_;
365 int ** Indices = Indices_;
366 double ** Values = Values_;
367 int NumMyCols_ = NumMyCols();
368
369 // If upper, point to last row
370 if ((Upper && !Trans) || (!Upper && Trans)) {
371 NumEntriesPerRow += NumMyRows_-1;
372 Indices += NumMyRows_-1;
373 Values += NumMyRows_-1;
374 }
375
376 double *xp = (double*)x.Values();
377 double *yp = (double*)y.Values();
378
379 if (!Trans) {
380
381 if (Upper) {
382
383 j0 = 1;
384 if (NoDiagonal()) j0--; // Include first term if no diagonal
385 for (i=NumMyRows_-1; i >=0; i--) {
386 int NumEntries = *NumEntriesPerRow--;
387 int * RowIndices = *Indices--;
388 double * RowValues = *Values--;
389 double sum = 0.0;
390 for (j=j0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
391
392 if (UnitDiagonal) yp[i] = xp[i] - sum;
393 else yp[i] = (xp[i] - sum)/RowValues[0];
394
395 }
396 }
397 else {
398 j0 = 1;
399 if (NoDiagonal()) j0--; // Include first term if no diagonal
400 for (i=0; i < NumMyRows_; i++) {
401 int NumEntries = *NumEntriesPerRow++ - j0;
402 int * RowIndices = *Indices++;
403 double * RowValues = *Values++;
404 double sum = 0.0;
405 for (j=0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
406
407 if (UnitDiagonal) yp[i] = xp[i] - sum;
408 else yp[i] = (xp[i] - sum)/RowValues[NumEntries];
409
410 }
411 }
412 }
413
414 // *********** Transpose case *******************************
415
416 else {
417
418 if (xp!=yp) for (i=0; i < NumMyCols_; i++) yp[i] = xp[i]; // Initialize y for transpose solve
419
420 if (Upper) {
421
422 j0 = 1;
423 if (NoDiagonal()) j0--; // Include first term if no diagonal
424
425 for (i=0; i < NumMyRows_; i++) {
426 int NumEntries = *NumEntriesPerRow++;
427 int * RowIndices = *Indices++;
428 double * RowValues = *Values++;
429 if (!UnitDiagonal) yp[i] = yp[i]/RowValues[0];
430 for (j=j0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
431 }
432 }
433 else {
434
435 j0 = 1;
436 if (NoDiagonal()) j0--; // Include first term if no diagonal
437
438 for (i=NumMyRows_-1; i >= 0; i--) {
439 int NumEntries = *NumEntriesPerRow-- - j0;
440 int * RowIndices = *Indices--;
441 double * RowValues = *Values--;
442 if (!UnitDiagonal) yp[i] = yp[i]/RowValues[NumEntries];
443 for (j=0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
444 }
445 }
446
447 }
448 UpdateFlops(2*NumGlobalNonzeros64());
449 return(0);
450}
451
452//=============================================================================
453int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
454 //
455 // This function find Y such that LY = X or UY = X or the transpose cases.
456 //
457 if (X.NumVectors()==1 && Y.NumVectors()==1) {
458 double * xp = (double *) X[0];
459 double * yp = (double *) Y[0];
460 Epetra_Vector x(View, X.Map(), xp);
461 Epetra_Vector y(View, Y.Map(), yp);
462 return(Solve(Upper, Trans, UnitDiagonal, x, y));
463 }
464 if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
465
466 if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
467 if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
468 if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
469 if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
470
471 int i, j, j0, k;
472 int * NumEntriesPerRow = NumEntriesPerRow_;
473 int ** Indices = Indices_;
474 double ** Values = Values_;
475 double diag;
476
477 // If upper, point to last row
478 if ((Upper && !Trans) || (!Upper && Trans)) {
479 NumEntriesPerRow += NumMyRows_-1;
480 Indices += NumMyRows_-1;
481 Values += NumMyRows_-1;
482 }
483
484 double **Xp = (double**)X.Pointers();
485 double **Yp = (double**)Y.Pointers();
486
487 int NumVectors = X.NumVectors();
488
489 if (!Trans) {
490
491
492 if (Upper) {
493
494 j0 = 1;
495 if (NoDiagonal()) j0--; // Include first term if no diagonal
496 for (i=NumMyRows_-1; i >=0; i--) {
497 int NumEntries = *NumEntriesPerRow--;
498 int * RowIndices = *Indices--;
499 double * RowValues = *Values--;
500 if (!UnitDiagonal) diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
501 for (k=0; k<NumVectors; k++) {
502 double sum = 0.0;
503 for (j=j0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
504
505 if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
506 else Yp[k][i] = (Xp[k][i] - sum)*diag;
507 }
508 }
509 }
510 else {
511 j0 = 1;
512 if (NoDiagonal()) j0--; // Include first term if no diagonal
513 for (i=0; i < NumMyRows_; i++) {
514 int NumEntries = *NumEntriesPerRow++ - j0;
515 int * RowIndices = *Indices++;
516 double * RowValues = *Values++;
517 if (!UnitDiagonal) diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
518 for (k=0; k<NumVectors; k++) {
519 double sum = 0.0;
520 for (j=0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
521
522 if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
523 else Yp[k][i] = (Xp[k][i] - sum)*diag;
524 }
525 }
526 }
527 }
528 // *********** Transpose case *******************************
529
530 else {
531
532 for (k=0; k<NumVectors; k++)
533 if (Yp[k]!=Xp[k]) for (i=0; i < NumMyRows_; i++) Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
534
535 if (Upper) {
536
537 j0 = 1;
538 if (NoDiagonal()) j0--; // Include first term if no diagonal
539
540 for (i=0; i < NumMyRows_; i++) {
541 int NumEntries = *NumEntriesPerRow++;
542 int * RowIndices = *Indices++;
543 double * RowValues = *Values++;
544 if (!UnitDiagonal) diag = 1.0/RowValues[j0]; // Take inverse of diagonal once for later use
545 for (k=0; k<NumVectors; k++) {
546 if (!UnitDiagonal) Yp[k][i] = Yp[k][i]*diag;
547 for (j=j0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
548 }
549 }
550 }
551 else {
552
553 j0 = 1;
554 if (NoDiagonal()) j0--; // Include first term if no diagonal
555
556 for (i=NumMyRows_-1; i>=0; i--) {
557 int NumEntries = *NumEntriesPerRow-- - j0;
558 int * RowIndices = *Indices--;
559 double * RowValues = *Values--;
560 for (k=0; k<NumVectors; k++) {
561 if (!UnitDiagonal) Yp[k][i] = Yp[k][i]/Xp[k][i];
562 for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
563 }
564 }
565 }
566 }
567
568 UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
569 return(0);
570}
@ Insert
#define EPETRA_CHK_ERR(a)
@ View
@ Copy
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
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_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.
double * Values() const
Get pointer to MultiVector values.
double ** Pointers() const
Get pointer to individual vector pointers.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.