Epetra Package Browser (Single Doxygen Collection)  Development
Epetra_JadMatrix.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 #include "Epetra_JadMatrix.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Import.h"
46 #include "Epetra_Export.h"
47 #include "Epetra_Vector.h"
48 #include "Epetra_MultiVector.h"
49 #include "Epetra_Comm.h"
50 #include "Epetra_Util.h"
51 #include "Epetra_IntSerialDenseVector.h"
52 #include "Epetra_RowMatrix.h"
53 #include "Epetra_CrsMatrix.h"
54 
55 //==============================================================================
56 Epetra_JadMatrix::Epetra_JadMatrix(const Epetra_RowMatrix & Matrix)
57  : Epetra_BasicRowMatrix(Matrix.RowMatrixRowMap().Comm()),
58  Values_(0),
59  Indices_(0),
60  IndexOffset_(0),
61  Profile_(0),
62  RowPerm_(0),
63  InvRowPerm_(0),
64  NumJaggedDiagonals_(Matrix.MaxNumEntries())
65 {
66  SetMaps(Matrix.RowMatrixRowMap(), Matrix.RowMatrixColMap(), Matrix.OperatorDomainMap(), Matrix.OperatorRangeMap());
67  if (!Matrix.Filled()) throw Matrix.RowMatrixRowMap().ReportError("Input matrix must have called FillComplete()", -1);
68  Allocate(Matrix);
69  SetLabel("Epetra::JadMatrix");
70 }
71 
72 //==============================================================================
73 Epetra_JadMatrix::~Epetra_JadMatrix(){}
74 
75 //==============================================================================
76 int Epetra_JadMatrix::UpdateValues(const Epetra_RowMatrix & Matrix, bool CheckStructure) {
77 
78  int NumEntries;
79  int * Indices = 0;
80  double * Values =0;
81 
82  int ierr = 0;
83 
84  try { // If matrix is an Epetra_CrsMatrix, we can get date much more cheaply
85 
86  const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Matrix);
87 
88  for (int i1=0; i1<NumMyRows_; i1++) {
89 
90  EPETRA_CHK_ERR(A.ExtractMyRowView(i1, NumEntries, Values, Indices)); // Get the current row based on the permutation
91  int i = InvRowPerm_[i1]; // Determine permuted row location
92  for (int j=0; j< NumEntries; j++) Values_[IndexOffset_[j]+i] = Values[j];
93  if (CheckStructure)
94  for (int j=0; j< NumEntries; j++) if (Indices_[IndexOffset_[j]+i] != Indices[j]) ierr = - 1;
95  }
96  }
97  catch (...) { // Otherwise just live with RowMatrix interface
98 
99  Epetra_SerialDenseVector curValues(NumJaggedDiagonals_);
100  Epetra_IntSerialDenseVector curIndices(NumJaggedDiagonals_);
101  Indices = curIndices.Values();
102  Values = curValues.Values();
103  for (int i1=0; i1<NumMyRows_; i1++) {
104  EPETRA_CHK_ERR(Matrix.ExtractMyRowCopy(i1, NumJaggedDiagonals_, NumEntries, Values, Indices)); // Get current row based on the permutation
105  int i = InvRowPerm_[i1]; // Determine permuted row location
106  for (int j=0; j< NumEntries; j++) Values_[IndexOffset_[j]+i] = Values[j];
107  if (CheckStructure)
108  for (int j=0; j< NumEntries; j++) if (Indices_[IndexOffset_[j]+i] != Indices[j]) ierr = - 1;
109  }
110  }
111 
112  HaveNumericConstants_ = false;
113  EPETRA_CHK_ERR(ierr);
114  return(ierr);
115 }
116 
117 //==============================================================================
118 void Epetra_JadMatrix::Allocate(const Epetra_RowMatrix & Matrix) {
119 
120  // Allocate IndexOffset storage
121  int numMyRows = Matrix.NumMyRows();
122  int numMyNonzeros = Matrix.NumMyNonzeros();
123 
124  IndexOffset_.Resize(NumJaggedDiagonals_+1);
125 
126  // Next compute permutation of rows
127  RowPerm_.Resize(numMyRows);
128  InvRowPerm_.Resize(numMyRows);
129  Profile_.Resize(numMyRows);
130  for (int i=0; i<numMyRows; i++) {
131  int NumEntries;
132  Matrix.NumMyRowEntries(i, NumEntries);
133  Profile_[i] = NumEntries;
134  RowPerm_[i] = i;
135  }
136 
137  Epetra_Util sorter;
138  int * RowPerm = RowPerm_.Values();
139  sorter.Sort(false, numMyRows, Profile_.Values(), 0, 0, 1, &RowPerm, 0, 0);
140  //cout << "Profile = " << Profile_ << std::endl;
141  //cout << "RowPerm = " << RowPerm_ << std::endl;
142  for (int i=0; i<numMyRows; i++) InvRowPerm_[RowPerm[i]] = i; // Compute inverse row permutation
143  //cout << "InvRowPerm = " << InvRowPerm_ << std::endl;
144 
145  // Now build IndexOffsets: These contain the lengths of the jagged diagonals
146 
147  for (int i=0; i<NumJaggedDiagonals_; i++) IndexOffset_[i] = 0;
148 
149  int curOffset = numMyRows;
150  int * curIndex = IndexOffset_.Values(); // point to first index (will be incremented immediately below)
151  for (int i=1; i<NumJaggedDiagonals_+1; i++) {
152  curIndex++;
153  while (*curIndex==0) {
154  if (Profile_[curOffset-1]<i) curOffset--;
155  else *curIndex = *(curIndex-1) + curOffset; // Set the length of the current jagged diagonal (plus scan sum)
156  }
157  }
158 
159  Values_.Resize(numMyNonzeros);
160  Indices_.Resize(numMyNonzeros);
161 
162  int NumEntries;
163  int * Indices = 0;
164  double * Values =0;
165 
166  try { // If matrix is an Epetra_CrsMatrix, we can get data much more cheaply
167 
168  const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Matrix);
169 
170  for (int i1=0; i1<numMyRows; i1++) {
171 
172  A.ExtractMyRowView(i1, NumEntries, Values, Indices); // Get the current row
173  int i = InvRowPerm_[i1]; // Determine permuted row location
174  //cout << "i1, i, NumEntries = " << i1 <<" "<< i <<" "<< NumEntries << std::endl;
175  for (int j=0; j< NumEntries; j++) {
176  Values_[IndexOffset_[j]+i] = Values[j];
177  Indices_[IndexOffset_[j]+i] = Indices[j];
178  }
179  }
180  }
181  catch (...) { // Otherwise just live with RowMatrix interface
182 
183  Epetra_SerialDenseVector curValues(NumJaggedDiagonals_);
184  Epetra_IntSerialDenseVector curIndices(NumJaggedDiagonals_);
185  Indices = curIndices.Values();
186  Values = curValues.Values();
187  for (int i1=0; i1<numMyRows; i1++) {
188  Matrix.ExtractMyRowCopy(i1, NumJaggedDiagonals_, NumEntries, Values, Indices); // Get current row based on the permutation
189  int i = InvRowPerm_[i1]; // Determine permuted row location
190  for (int j=0; j< NumEntries; j++) {
191  Values_[IndexOffset_[j]+i] = Values[j];
192  Indices_[IndexOffset_[j]+i] = Indices[j];
193  }
194  }
195  }
196 }
197 //=============================================================================
198 int Epetra_JadMatrix::NumMyRowEntries(int MyRow, int & NumEntries) const {
199  int i = InvRowPerm_[MyRow]; // Determine permuted row location
200  NumEntries = Profile_[i]; // NNZ in current row
201  return(0);
202 }
203 //=============================================================================
204 int Epetra_JadMatrix::ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const {
205 
206  if(MyRow < 0 || MyRow >= NumMyRows_)
207  EPETRA_CHK_ERR(-1); // Not in Row range
208 
209  int i = InvRowPerm_[MyRow]; // Determine permuted row location
210  NumEntries = Profile_[i]; // NNZ in current row
211  if(NumEntries > Length)
212  EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumEntries
213 
214  for (int j=0; j< NumEntries; j++) Values[j] = Values_[IndexOffset_[j]+i];
215  for (int j=0; j< NumEntries; j++) Indices[j] = Indices_[IndexOffset_[j]+i];
216  return(0);
217 }
218 //=============================================================================
219 int Epetra_JadMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
220  //
221  // This function forms the product Y = A * Y or Y = A' * X
222  //
223 
224  int NumVectors = X.NumVectors();
225  if (NumVectors!=Y.NumVectors()) {
226  EPETRA_CHK_ERR(-1); // Need same number of vectors in each MV
227  }
228 
229  double** Xp = (double**) X.Pointers();
230  double** Yp = (double**) Y.Pointers();
231  int LDX = X.ConstantStride() ? X.Stride() : 0;
232  int LDY = Y.ConstantStride() ? Y.Stride() : 0;
233  UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
234  UpdateExportVector(NumVectors);
235 
236  if (!TransA) {
237 
238  // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
239  if (Importer()!=0) {
240  EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
241  Xp = (double**)ImportVector_->Pointers();
242  LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
243  }
244 
245  // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
246  if (Exporter()!=0) {
247  Yp = (double**)ExportVector_->Pointers();
248  LDY = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
249  }
250 
251  // Do actual computation
252  if (NumVectors==1)
253  GeneralMV(TransA, *Xp, *Yp);
254  else
255  GeneralMM(TransA, Xp, LDX, Yp, LDY, NumVectors);
256  if (Exporter()!=0) {
257  Y.PutScalar(0.0); // Make sure target is zero
258  Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
259  }
260  // Handle case of rangemap being a local replicated map
261  if (!OperatorRangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
262  }
263  else { // Transpose operation
264 
265 
266  // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
267 
268  if (Exporter()!=0) {
269  EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
270  Xp = (double**)ExportVector_->Pointers();
271  LDX = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
272  }
273 
274  // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
275  if (Importer()!=0) {
276  Yp = (double**)ImportVector_->Pointers();
277  LDY = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
278  }
279 
280  // Do actual computation
281  if (NumVectors==1)
282  GeneralMV(TransA, *Xp, *Yp);
283  else
284  GeneralMM(TransA, Xp, LDX, Yp, LDY, NumVectors);
285  if (Importer()!=0) {
286  Y.PutScalar(0.0); // Make sure target is zero
287  EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
288  }
289  // Handle case of rangemap being a local replicated map
290  if (!OperatorDomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
291  }
292 
293  UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
294  return(0);
295 }
296 //=======================================================================================================
297 void Epetra_JadMatrix::GeneralMM(bool TransA, double ** X, int LDX, double ** Y, int LDY, int NumVectors) const {
298 
299  if (LDX==0 || LDY==0 || NumVectors==1) {// Can't unroll RHS if X or Y not strided
300  for (int k=0; k<NumVectors; k++) GeneralMV(TransA, X[k], Y[k]);
301  }
302  else if (NumVectors==2) // Special 2 RHS case (does unrolling in both NumVectors and NumJaggedDiagonals)
303  GeneralMM2RHS(TransA, X[0], LDX, Y[0], LDY);
304  // Otherwise unroll RHS only
305  else
306  GeneralMM3RHS(TransA, X, LDX, Y, LDY, NumVectors);
307 
308  return;
309 }
310 //=======================================================================================================
311 void Epetra_JadMatrix::GeneralMM3RHS(bool TransA, double ** X, int ldx, double ** Y, int ldy, int NumVectors) const {
312 
313 #ifdef _CRAY
314 #define Pragma(S) _Pragma(S)
315 #else
316 #define Pragma(S)
317 #endif
318 
319  // Routine for 3 or more RHS
320 
321  const double * Values = Values_.Values();
322  const int * Indices = Indices_.Values();
323  const int * IndexOffset = IndexOffset_.Values();
324  const int * RowPerm = RowPerm_.Values();
325  for (int j=0; j<NumVectors; j++) {
326  double * y = Y[j];
327  if (!TransA)
328  for (int i=0; i<NumMyRows_; i++) y[i] = 0.0;
329  else
330  for (int i=0; i<NumMyCols_; i++) y[i] = 0.0;
331  }
332 
333  int nv = NumVectors%5; if (nv==0) nv=5;
334  double * x = X[0];
335  double * y = Y[0];
336 
337 
338  for (int k=0; k<NumVectors; k+=5) {
339 
340  for (int j=0; j<NumJaggedDiagonals_; j++) {
341  const int * curIndices = Indices+IndexOffset[j];
342  const double * curValues = Values+IndexOffset[j];
343  int jaggedDiagonalLength = IndexOffset[j+1]-IndexOffset[j];
344  switch (nv){
345  case 1:
346  {
347  if (!TransA) {
348 Pragma("_CRI ivdep")
349  for (int i=0; i<jaggedDiagonalLength; i++) {
350  int ix = curIndices[i];
351  int iy = RowPerm[i];
352  double val = curValues[i];
353  y[iy] += val*x[ix];
354  }
355  }
356  else {
357 Pragma("_CRI ivdep")
358  for (int i=0; i<jaggedDiagonalLength; i++) {
359  int iy = curIndices[i];
360  int ix = RowPerm[i];
361  double val = curValues[i];
362  y[iy] += val*x[ix];
363  }
364  }
365  break;
366  }
367  case 2:
368  {
369  if (!TransA) {
370 Pragma("_CRI ivdep")
371  for (int i=0; i<jaggedDiagonalLength; i++) {
372  int ix = curIndices[i];
373  int iy = RowPerm[i];
374  double val = curValues[i];
375  y[iy] += val*x[ix];
376  iy+=ldy; ix+=ldx;
377  y[iy] += val*x[ix];
378  }
379  }
380  else {
381 Pragma("_CRI ivdep")
382  for (int i=0; i<jaggedDiagonalLength; i++) {
383  int iy = curIndices[i];
384  int ix = RowPerm[i];
385  double val = curValues[i];
386  y[iy] += val*x[ix];
387  iy+=ldy; ix+=ldx;
388  y[iy] += val*x[ix];
389  }
390  }
391  break;
392  }
393  case 3:
394  {
395  if (!TransA) {
396 Pragma("_CRI ivdep")
397  for (int i=0; i<jaggedDiagonalLength; i++) {
398  int ix = curIndices[i];
399  int iy = RowPerm[i];
400  double val = curValues[i];
401  y[iy] += val*x[ix];
402  iy+=ldy; ix+=ldx;
403  y[iy] += val*x[ix];
404  iy+=ldy; ix+=ldx;
405  y[iy] += val*x[ix];
406  }
407  }
408  else {
409 Pragma("_CRI ivdep")
410  for (int i=0; i<jaggedDiagonalLength; i++) {
411  int iy = curIndices[i];
412  int ix = RowPerm[i];
413  double val = curValues[i];
414  y[iy] += val*x[ix];
415  iy+=ldy; ix+=ldx;
416  y[iy] += val*x[ix];
417  iy+=ldy; ix+=ldx;
418  y[iy] += val*x[ix];
419  }
420  }
421  break;
422  }
423  case 4:
424  {
425  if (!TransA) {
426 Pragma("_CRI ivdep")
427  for (int i=0; i<jaggedDiagonalLength; i++) {
428  int ix = curIndices[i];
429  int iy = RowPerm[i];
430  double val = curValues[i];
431  y[iy] += val*x[ix];
432  iy+=ldy; ix+=ldx;
433  y[iy] += val*x[ix];
434  iy+=ldy; ix+=ldx;
435  y[iy] += val*x[ix];
436  iy+=ldy; ix+=ldx;
437  y[iy] += val*x[ix];
438  }
439  }
440  else {
441 Pragma("_CRI ivdep")
442  for (int i=0; i<jaggedDiagonalLength; i++) {
443  int iy = curIndices[i];
444  int ix = RowPerm[i];
445  double val = curValues[i];
446  y[iy] += val*x[ix];
447  iy+=ldy; ix+=ldx;
448  y[iy] += val*x[ix];
449  iy+=ldy; ix+=ldx;
450  y[iy] += val*x[ix];
451  iy+=ldy; ix+=ldx;
452  y[iy] += val*x[ix];
453  }
454  }
455  break;
456  }
457  case 5:
458  {
459  if (!TransA) {
460 Pragma("_CRI ivdep")
461  for (int i=0; i<jaggedDiagonalLength; i++) {
462  int ix = curIndices[i];
463  int iy = RowPerm[i];
464  double val = curValues[i];
465  y[iy] += val*x[ix];
466  iy+=ldy; ix+=ldx;
467  y[iy] += val*x[ix];
468  iy+=ldy; ix+=ldx;
469  y[iy] += val*x[ix];
470  iy+=ldy; ix+=ldx;
471  y[iy] += val*x[ix];
472  iy+=ldy; ix+=ldx;
473  y[iy] += val*x[ix];
474  }
475  }
476  else {
477 Pragma("_CRI ivdep")
478  for (int i=0; i<jaggedDiagonalLength; i++) {
479  int iy = curIndices[i];
480  int ix = RowPerm[i];
481  double val = curValues[i];
482  y[iy] += val*x[ix];
483  iy+=ldy; ix+=ldx;
484  y[iy] += val*x[ix];
485  iy+=ldy; ix+=ldx;
486  y[iy] += val*x[ix];
487  iy+=ldy; ix+=ldx;
488  y[iy] += val*x[ix];
489  iy+=ldy; ix+=ldx;
490  y[iy] += val*x[ix];
491  }
492  }
493  break;
494  }
495  }
496  }
497  x += nv*ldx;
498  y += nv*ldy;
499  nv = 5; // After initial remainder, we will always do 5 RHS
500  }
501  return;
502 }
503 //=======================================================================================================
504 void Epetra_JadMatrix::GeneralMM2RHS(bool TransA, double * x, int ldx, double * y, int ldy) const {
505 
506  // special 2 rhs case
507 
508  const double * Values = Values_.Values();
509  const int * Indices = Indices_.Values();
510  const int * IndexOffset = IndexOffset_.Values();
511  const int * RowPerm = RowPerm_.Values();
512  if (!TransA)
513  for (int i=0; i<NumMyRows_; i++) {
514  y[i] = 0.0;
515  y[i+ldy] = 0.0;
516  }
517  else
518  for (int i=0; i<NumMyCols_; i++) {
519  y[i] = 0.0;
520  y[i+ldy] = 0.0;
521  }
522 
523  int j = 0;
524  while (j<NumJaggedDiagonals_) {
525  int j0 = j;
526  int jaggedDiagonalLength = IndexOffset[j+1]-IndexOffset[j];
527  j++;
528  // check if other diagonals have same length up to a max of 2
529  while ((j<NumJaggedDiagonals_-1) && (IndexOffset[j+1]-IndexOffset[j]==jaggedDiagonalLength) && (j-j0<2)) j++;
530 
531  int numDiags = j-j0;
532  assert(numDiags<3 && numDiags>0);
533  assert(j<NumJaggedDiagonals_+1);
534 
535  switch (numDiags){
536  case 1:
537  {
538  const int * curIndices = Indices+IndexOffset[j0];
539  const double * curValues = Values+IndexOffset[j0];
540  if (!TransA) {
541 Pragma("_CRI ivdep")
542  for (int i=0; i<jaggedDiagonalLength; i++) {
543  int ix = curIndices[i];
544  int iy = RowPerm[i];
545  y[iy] += curValues[i]*x[ix];
546  iy+=ldy; ix+=ldx;
547  y[iy] += curValues[i]*x[ix];
548  }
549  }
550  else {
551 Pragma("_CRI ivdep")
552  for (int i=0; i<jaggedDiagonalLength; i++){
553  int iy = curIndices[i];
554  int ix = RowPerm[i];
555  y[iy] += curValues[i]*x[ix];
556  iy+=ldy; ix+=ldx;
557  y[iy] += curValues[i]*x[ix];
558  }
559  }
560  break;
561  }
562  case 2:
563  {
564  const int * curIndices0 = Indices+IndexOffset[j0];
565  const double * curValues0 = Values+IndexOffset[j0++];
566  const int * curIndices1 = Indices+IndexOffset[j0];
567  const double * curValues1 = Values+IndexOffset[j0];
568  if (!TransA) {
569 Pragma("_CRI ivdep")
570  for (int i=0; i<jaggedDiagonalLength; i++) {
571  int ix0 = curIndices0[i];
572  int ix1 = curIndices1[i];
573  int iy = RowPerm[i];
574  y[iy] +=
575  curValues0[i]*x[ix0] +
576  curValues1[i]*x[ix1];
577  iy+=ldy; ix0+=ldx; ix1+=ldx;
578  y[iy] +=
579  curValues0[i]*x[ix0] +
580  curValues1[i]*x[ix1];
581  }
582  }
583  else {
584 Pragma("_CRI ivdep")
585  for (int i=0; i<jaggedDiagonalLength; i++) {
586  int iy0 = curIndices0[i];
587  int iy1 = curIndices1[i];
588  int ix = RowPerm[i];
589  double xval = x[ix];
590  y[iy0] += curValues0[i]*xval;
591  y[iy1] += curValues1[i]*xval;
592  ix+=ldx; iy0+=ldy; iy1+=ldy;
593  xval = x[ix];
594  y[iy0] += curValues0[i]*xval;
595  y[iy1] += curValues1[i]*xval;
596  }
597  }
598  }
599  break;
600  }
601  }
602  return;
603 }
604 //=======================================================================================================
605 void Epetra_JadMatrix::GeneralMV(bool TransA, double * x, double * y) const {
606 
607  const double * Values = Values_.Values();
608  const int * Indices = Indices_.Values();
609  const int * IndexOffset = IndexOffset_.Values();
610  const int * RowPerm = RowPerm_.Values();
611  if (!TransA)
612  for (int i=0; i<NumMyRows_; i++) y[i] = 0.0;
613  else
614  for (int i=0; i<NumMyCols_; i++) y[i] = 0.0;
615 
616  int j = 0;
617  while (j<NumJaggedDiagonals_) {
618  int j0 = j;
619  int jaggedDiagonalLength = IndexOffset[j+1]-IndexOffset[j];
620  j++;
621  // check if other diagonals have same length up to a max of 5
622  while ((j<NumJaggedDiagonals_-1) && (IndexOffset[j+1]-IndexOffset[j]==jaggedDiagonalLength) && (j-j0<5)) j++;
623 
624  int numDiags = j-j0;
625  assert(numDiags<6 && numDiags>0);
626  assert(j<NumJaggedDiagonals_+1);
627 
628  switch (numDiags){
629  case 1:
630  {
631  const int * curIndices = Indices+IndexOffset[j0];
632  const double * curValues = Values+IndexOffset[j0];
633  if (!TransA) {
634 Pragma("_CRI ivdep")
635  for (int i=0; i<jaggedDiagonalLength; i++)
636  y[RowPerm[i]] += curValues[i]*x[curIndices[i]];
637  }
638  else {
639 Pragma("_CRI ivdep")
640  for (int i=0; i<jaggedDiagonalLength; i++)
641  y[curIndices[i]] += curValues[i]*x[RowPerm[i]];
642  }
643  break;
644  }
645  case 2:
646  {
647  const int * curIndices0 = Indices+IndexOffset[j0];
648  const double * curValues0 = Values+IndexOffset[j0++];
649  const int * curIndices1 = Indices+IndexOffset[j0];
650  const double * curValues1 = Values+IndexOffset[j0];
651  if (!TransA) {
652 Pragma("_CRI ivdep")
653  for (int i=0; i<jaggedDiagonalLength; i++) {
654  y[RowPerm[i]] +=
655  curValues0[i]*x[curIndices0[i]] +
656  curValues1[i]*x[curIndices1[i]];
657  }
658  }
659  else {
660  //Pragma("_CRI ivdep") (Collisions possible)
661  for (int i=0; i<jaggedDiagonalLength; i++) {
662  double xval = x[RowPerm[i]];
663  y[curIndices0[i]] += curValues0[i]*xval;
664  y[curIndices1[i]] += curValues1[i]*xval;
665  }
666  }
667  }
668  break;
669  case 3:
670  {
671  const int * curIndices0 = Indices+IndexOffset[j0];
672  const double * curValues0 = Values+IndexOffset[j0++];
673  const int * curIndices1 = Indices+IndexOffset[j0];
674  const double * curValues1 = Values+IndexOffset[j0++];
675  const int * curIndices2 = Indices+IndexOffset[j0];
676  const double * curValues2 = Values+IndexOffset[j0];
677  if (!TransA) {
678 Pragma("_CRI ivdep")
679  for (int i=0; i<jaggedDiagonalLength; i++) {
680  y[RowPerm[i]] +=
681  curValues0[i]*x[curIndices0[i]] +
682  curValues1[i]*x[curIndices1[i]] +
683  curValues2[i]*x[curIndices2[i]];
684  }
685  }
686  else {
687  //Pragma("_CRI ivdep") (Collisions possible)
688  for (int i=0; i<jaggedDiagonalLength; i++) {
689  double xval = x[RowPerm[i]];
690  y[curIndices0[i]] += curValues0[i]*xval;
691  y[curIndices1[i]] += curValues1[i]*xval;
692  y[curIndices2[i]] += curValues2[i]*xval;
693  }
694  }
695  }
696  break;
697  case 4:
698  {
699  const int * curIndices0 = Indices+IndexOffset[j0];
700  const double * curValues0 = Values+IndexOffset[j0++];
701  const int * curIndices1 = Indices+IndexOffset[j0];
702  const double * curValues1 = Values+IndexOffset[j0++];
703  const int * curIndices2 = Indices+IndexOffset[j0];
704  const double * curValues2 = Values+IndexOffset[j0++];
705  const int * curIndices3 = Indices+IndexOffset[j0];
706  const double * curValues3 = Values+IndexOffset[j0];
707  if (!TransA) {
708 Pragma("_CRI ivdep")
709  for (int i=0; i<jaggedDiagonalLength; i++) {
710  y[RowPerm[i]] +=
711  curValues0[i]*x[curIndices0[i]] +
712  curValues1[i]*x[curIndices1[i]] +
713  curValues2[i]*x[curIndices2[i]] +
714  curValues3[i]*x[curIndices3[i]];
715  }
716  }
717  else {
718  //Pragma("_CRI ivdep") (Collisions possible)
719  for (int i=0; i<jaggedDiagonalLength; i++) {
720  double xval = x[RowPerm[i]];
721  y[curIndices0[i]] += curValues0[i]*xval;
722  y[curIndices1[i]] += curValues1[i]*xval;
723  y[curIndices2[i]] += curValues2[i]*xval;
724  y[curIndices3[i]] += curValues3[i]*xval;
725  }
726  }
727  }
728  break;
729  case 5:
730  {
731  const int * curIndices0 = Indices+IndexOffset[j0];
732  const double * curValues0 = Values+IndexOffset[j0++];
733  const int * curIndices1 = Indices+IndexOffset[j0];
734  const double * curValues1 = Values+IndexOffset[j0++];
735  const int * curIndices2 = Indices+IndexOffset[j0];
736  const double * curValues2 = Values+IndexOffset[j0++];
737  const int * curIndices3 = Indices+IndexOffset[j0];
738  const double * curValues3 = Values+IndexOffset[j0++];
739  const int * curIndices4 = Indices+IndexOffset[j0];
740  const double * curValues4 = Values+IndexOffset[j0];
741  if (!TransA) {
742 Pragma("_CRI ivdep")
743  for (int i=0; i<jaggedDiagonalLength; i++) {
744  y[RowPerm[i]] +=
745  curValues0[i]*x[curIndices0[i]] +
746  curValues1[i]*x[curIndices1[i]] +
747  curValues2[i]*x[curIndices2[i]] +
748  curValues3[i]*x[curIndices3[i]] +
749  curValues4[i]*x[curIndices4[i]];
750  }
751  }
752  else {
753  // Pragma("_CRI ivdep") (Collisions possible)
754  for (int i=0; i<jaggedDiagonalLength; i++) {
755  double xval = x[RowPerm[i]];
756  y[curIndices0[i]] += curValues0[i]*xval;
757  y[curIndices1[i]] += curValues1[i]*xval;
758  y[curIndices2[i]] += curValues2[i]*xval;
759  y[curIndices3[i]] += curValues3[i]*xval;
760  y[curIndices4[i]] += curValues4[i]*xval;
761  }
762  }
763  }
764  break;
765  }
766  }
767  return;
768 }
Epetra_BasicRowMatrix
Epetra_BasicRowMatrix: A class for simplifying the development of Epetra_RowMatrix adapters.
Definition: Epetra_BasicRowMatrix.h:110
Epetra_SerialDenseVector
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
Definition: Epetra_SerialDenseVector.h:95
Epetra_Util::Sort
static void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_Util Sort Routine (Shell sort)
Definition: Epetra_Util.cpp:114
Epetra_IntSerialDenseVector.h
Epetra_Comm::NumProc
virtual int NumProc() const =0
Returns total number of processes.
Epetra_JadMatrix::Indices_
Epetra_IntSerialDenseVector Indices_
Definition: Epetra_JadMatrix.h:221
Add
Definition: Epetra_CombineMode.h:64
Epetra_MultiVector::PutScalar
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Definition: Epetra_MultiVector.cpp:595
Epetra_JadMatrix::ExtractMyRowCopy
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.
Definition: Epetra_JadMatrix.cpp:204
Epetra_CompObject::UpdateFlops
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Definition: Epetra_CompObject.h:99
EPETRA_CHK_ERR
#define EPETRA_CHK_ERR(a)
Definition: Epetra_ConfigDefs.h:307
Epetra_RowMatrix.h
Epetra_RowMatrix::NumMyRowEntries
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
Returns the number of nonzero entries in MyRow.
Epetra_BasicRowMatrix::NumMyCols_
int NumMyCols_
Definition: Epetra_BasicRowMatrix.h:548
Epetra_CrsMatrix.h
Epetra_JadMatrix::GeneralMM3RHS
void GeneralMM3RHS(bool TransA, double **X, int LDX, double **Y, int LDY, int NumVectors) const
Definition: Epetra_JadMatrix.cpp:311
Epetra_JadMatrix::Profile_
Epetra_IntSerialDenseVector Profile_
Definition: Epetra_JadMatrix.h:223
Epetra_JadMatrix.h
Epetra_IntSerialDenseVector::Values
int * Values()
Returns pointer to the values in vector.
Definition: Epetra_IntSerialDenseVector.h:211
Epetra_Vector.h
Epetra_BasicRowMatrix::NumMyRows_
int NumMyRows_
Definition: Epetra_BasicRowMatrix.h:547
Epetra_JadMatrix::GeneralMV
void GeneralMV(bool TransA, double *x, double *y) const
Definition: Epetra_JadMatrix.cpp:605
Epetra_RowMatrix::Filled
virtual bool Filled() const =0
If FillComplete() has been called, this query returns true, otherwise it returns false.
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
Epetra_MultiVector::ConstantStride
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
Definition: Epetra_MultiVector.h:944
Insert
Definition: Epetra_CombineMode.h:68
Epetra_JadMatrix::Multiply
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_JadMatrix multiplied by a Epetra_MultiVector X in Y.
Definition: Epetra_JadMatrix.cpp:219
Epetra_CrsMatrix
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
Definition: Epetra_CrsMatrix.h:173
Epetra_JadMatrix::IndexOffset_
Epetra_IntSerialDenseVector IndexOffset_
Definition: Epetra_JadMatrix.h:222
Epetra_JadMatrix::RowPerm_
Epetra_IntSerialDenseVector RowPerm_
Definition: Epetra_JadMatrix.h:224
Epetra_BasicRowMatrix::ExportVector_
Epetra_MultiVector * ExportVector_
Definition: Epetra_BasicRowMatrix.h:560
Epetra_RowMatrix::NumMyRows
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
Epetra_DistObject::Export
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.
Definition: Epetra_DistObject.cpp:198
Epetra_Util.h
Epetra_MultiVector::Pointers
double ** Pointers() const
Get pointer to individual vector pointers.
Definition: Epetra_MultiVector.h:988
Epetra_SerialDenseVector::Resize
int Resize(int Length_in)
Resize a Epetra_SerialDenseVector object.
Definition: Epetra_SerialDenseVector.h:171
Epetra_JadMatrix::~Epetra_JadMatrix
virtual ~Epetra_JadMatrix()
Epetra_JadMatrix Destructor.
Definition: Epetra_JadMatrix.cpp:73
Epetra_BasicRowMatrix::HaveNumericConstants_
bool HaveNumericConstants_
Definition: Epetra_BasicRowMatrix.h:555
Epetra_JadMatrix::NumMyRowEntries
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
Definition: Epetra_JadMatrix.cpp:198
Epetra_BasicRowMatrix::Importer
virtual const Epetra_Import * Importer() const
Returns the Epetra_Import object that contains the import operations for distributed operations,...
Definition: Epetra_BasicRowMatrix.h:499
Epetra_JadMatrix::InvRowPerm_
Epetra_IntSerialDenseVector InvRowPerm_
Definition: Epetra_JadMatrix.h:225
Epetra_RowMatrix
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
Definition: Epetra_RowMatrix.h:68
Epetra_JadMatrix::Allocate
void Allocate(const Epetra_RowMatrix &Matrix)
Definition: Epetra_JadMatrix.cpp:118
Epetra_JadMatrix::Epetra_JadMatrix
Epetra_JadMatrix(const Epetra_RowMatrix &Matrix)
Epetra_JadMatrix constuctor.
Definition: Epetra_JadMatrix.cpp:56
Epetra_Object::SetLabel
virtual void SetLabel(const char *const Label)
Epetra_Object Label definition using char *.
Definition: Epetra_Object.cpp:134
Epetra_MultiVector::Reduce
int Reduce()
Definition: Epetra_MultiVector.cpp:2400
Epetra_Object::ReportError
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
Definition: Epetra_Object.cpp:103
Epetra_Util
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
Epetra_Export.h
Epetra_BasicRowMatrix::OperatorRangeMap
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator (same as domain).
Definition: Epetra_BasicRowMatrix.h:409
Epetra_JadMatrix::NumJaggedDiagonals_
int NumJaggedDiagonals_
Definition: Epetra_JadMatrix.h:226
Epetra_IntSerialDenseVector::Resize
int Resize(int Length_in)
Resize a Epetra_IntSerialDenseVector object.
Definition: Epetra_IntSerialDenseVector.h:153
Epetra_JadMatrix::GeneralMM
void GeneralMM(bool TransA, double **X, int LDX, double **Y, int LDY, int NumVectors) const
Definition: Epetra_JadMatrix.cpp:297
Epetra_Import.h
Epetra_BlockMap::DistributedGlobal
bool DistributedGlobal() const
Returns true if map is defined across more than one processor.
Definition: Epetra_BlockMap.h:694
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
Returns a copy of the specified local row in user-provided arrays.
Epetra_JadMatrix::Values_
Epetra_SerialDenseVector Values_
Definition: Epetra_JadMatrix.h:220
Epetra_Comm.h
Epetra_RowMatrix::NumMyNonzeros
virtual int NumMyNonzeros() const =0
Returns the number of nonzero entries in the calling processor's portion of the matrix.
Epetra_IntSerialDenseVector
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
Definition: Epetra_IntSerialDenseVector.h:87
Epetra_BasicRowMatrix::OperatorDomainMap
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Definition: Epetra_BasicRowMatrix.h:406
Epetra_RowMatrix::RowMatrixColMap
virtual const Epetra_Map & RowMatrixColMap() const =0
Returns the Epetra_Map object associated with the columns of this matrix.
Epetra_BasicRowMatrix::ImportVector_
Epetra_MultiVector * ImportVector_
Definition: Epetra_BasicRowMatrix.h:559
Epetra_MultiVector
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
Definition: Epetra_MultiVector.h:184
A
Epetra_Operator::OperatorRangeMap
virtual const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
Epetra_BasicRowMatrix::SetMaps
void SetMaps(const Epetra_Map &RowMap, const Epetra_Map &ColMap)
Set maps (Version 1); call this function or the next, but not both.
Definition: Epetra_BasicRowMatrix.cpp:101
Epetra_Map.h
Epetra_JadMatrix::UpdateValues
int UpdateValues(const Epetra_RowMatrix &Matrix, bool CheckStructure=false)
Update values using a matrix with identical structure.
Definition: Epetra_JadMatrix.cpp:76
Epetra_DistObject::Import
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.
Definition: Epetra_DistObject.cpp:113
Epetra_BasicRowMatrix::UpdateImportVector
void UpdateImportVector(int NumVectors) const
Definition: Epetra_BasicRowMatrix.cpp:456
Epetra_BasicRowMatrix::Exporter
virtual const Epetra_Export * Exporter() const
Returns the Epetra_Export object that contains the export operations for distributed operations,...
Definition: Epetra_BasicRowMatrix.h:511
Epetra_Operator::OperatorDomainMap
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
Epetra_BasicRowMatrix::Comm
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
Definition: Epetra_BasicRowMatrix.h:424
Pragma
#define Pragma(S)
Epetra_BasicRowMatrix::UpdateExportVector
void UpdateExportVector(int NumVectors) const
Definition: Epetra_BasicRowMatrix.cpp:470
Epetra_MultiVector.h
Epetra_SerialDenseVector::Values
double * Values() const
Returns pointer to the values in vector.
Definition: Epetra_SerialDenseVector.h:274
Epetra_MultiVector::Stride
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true).
Definition: Epetra_MultiVector.h:941
Epetra_BasicRowMatrix::NumGlobalNonzeros64
virtual long long NumGlobalNonzeros64() const
Definition: Epetra_BasicRowMatrix.h:352
Epetra_JadMatrix::GeneralMM2RHS
void GeneralMM2RHS(bool TransA, double *x, int ldx, double *y, int ldy) const
Definition: Epetra_JadMatrix.cpp:504
Epetra_MultiVector::NumVectors
int NumVectors() const
Returns the number of vectors in the multi-vector.
Definition: Epetra_MultiVector.h:925
Generated by   doxygen 1.8.16