1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // Implementation of sparse matrix
13
14 #include "linalg.hpp"
15 #include "../general/forall.hpp"
16 #include "../general/table.hpp"
17 #include "../general/sort_pairs.hpp"
18
19 #include <iostream>
20 #include <iomanip>
21 #include <cmath>
22 #include <algorithm>
23 #include <limits>
24 #include <cstring>
25
26 namespace mfem
27 {
28
29 using namespace std;
30
31 #ifdef MFEM_USE_CUDA
32 int SparseMatrix::SparseMatrixCount = 0;
33 cusparseHandle_t SparseMatrix::handle = nullptr;
34 size_t SparseMatrix::bufferSize = 0;
35 void * SparseMatrix::dBuffer = nullptr;
36 #endif
37
InitCuSparse()38 void SparseMatrix::InitCuSparse()
39 {
40 // Initialize cuSPARSE library
41 #ifdef MFEM_USE_CUDA
42 if (Device::Allows(Backend::CUDA_MASK))
43 {
44 if (!handle) { cusparseCreate(&handle); }
45 useCuSparse=true;
46 SparseMatrixCount++;
47 }
48 else
49 {
50 useCuSparse=false;
51 }
52 #endif
53 }
54
ClearCuSparse()55 void SparseMatrix::ClearCuSparse()
56 {
57 #ifdef MFEM_USE_CUDA
58 if (initBuffers)
59 {
60 #if CUDA_VERSION >= 10010
61 cusparseDestroySpMat(matA_descr);
62 cusparseDestroyDnVec(vecX_descr);
63 cusparseDestroyDnVec(vecY_descr);
64 #else
65 cusparseDestroyMatDescr(matA_descr);
66 #endif
67 initBuffers = false;
68 }
69 #endif
70 }
71
SparseMatrix(int nrows,int ncols)72 SparseMatrix::SparseMatrix(int nrows, int ncols)
73 : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
74 Rows(new RowNode *[nrows]),
75 current_row(-1),
76 ColPtrJ(NULL),
77 ColPtrNode(NULL),
78 At(NULL),
79 isSorted(false)
80 {
81 // We probably do not need to set the ownership flags here.
82 I.Reset(); I.SetHostPtrOwner(true);
83 J.Reset(); J.SetHostPtrOwner(true);
84 A.Reset(); A.SetHostPtrOwner(true);
85
86 for (int i = 0; i < nrows; i++)
87 {
88 Rows[i] = NULL;
89 }
90
91 #ifdef MFEM_USE_MEMALLOC
92 NodesMem = new RowNodeAlloc;
93 #endif
94
95 InitCuSparse();
96 }
97
SparseMatrix(int * i,int * j,double * data,int m,int n)98 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
99 : AbstractSparseMatrix(m, n),
100 Rows(NULL),
101 ColPtrJ(NULL),
102 ColPtrNode(NULL),
103 At(NULL),
104 isSorted(false)
105 {
106 I.Wrap(i, height+1, true);
107 J.Wrap(j, I[height], true);
108 A.Wrap(data, I[height], true);
109
110 #ifdef MFEM_USE_MEMALLOC
111 NodesMem = NULL;
112 #endif
113
114 InitCuSparse();
115 }
116
SparseMatrix(int * i,int * j,double * data,int m,int n,bool ownij,bool owna,bool issorted)117 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
118 bool ownij, bool owna, bool issorted)
119 : AbstractSparseMatrix(m, n),
120 Rows(NULL),
121 ColPtrJ(NULL),
122 ColPtrNode(NULL),
123 At(NULL),
124 isSorted(issorted)
125 {
126 I.Wrap(i, height+1, ownij);
127 J.Wrap(j, I[height], ownij);
128
129 #ifdef MFEM_USE_MEMALLOC
130 NodesMem = NULL;
131 #endif
132
133 if (data)
134 {
135 A.Wrap(data, I[height], owna);
136 }
137 else
138 {
139 const int nnz = I[height];
140 A.New(nnz);
141 for (int i=0; i<nnz; ++i)
142 {
143 A[i] = 0.0;
144 }
145 }
146
147 InitCuSparse();
148 }
149
SparseMatrix(int nrows,int ncols,int rowsize)150 SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
151 : AbstractSparseMatrix(nrows, ncols)
152 , Rows(NULL)
153 , ColPtrJ(NULL)
154 , ColPtrNode(NULL)
155 , At(NULL)
156 , isSorted(false)
157 {
158 #ifdef MFEM_USE_MEMALLOC
159 NodesMem = NULL;
160 #endif
161 I.New(nrows + 1);
162 J.New(nrows * rowsize);
163 A.New(nrows * rowsize);
164
165 for (int i = 0; i <= nrows; i++)
166 {
167 I[i] = i * rowsize;
168 }
169
170 InitCuSparse();
171 }
172
SparseMatrix(const SparseMatrix & mat,bool copy_graph,MemoryType mt)173 SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph,
174 MemoryType mt)
175 : AbstractSparseMatrix(mat.Height(), mat.Width())
176 {
177 if (mat.Finalized())
178 {
179 const int nnz = mat.I[height];
180 if (copy_graph)
181 {
182 I.New(height+1, mt == MemoryType::PRESERVE ? mat.I.GetMemoryType() : mt);
183 J.New(nnz, mt == MemoryType::PRESERVE ? mat.J.GetMemoryType() : mt);
184 I.CopyFrom(mat.I, height+1);
185 J.CopyFrom(mat.J, nnz);
186 }
187 else
188 {
189 I = mat.I;
190 J = mat.J;
191 I.ClearOwnerFlags();
192 J.ClearOwnerFlags();
193 }
194 A.New(nnz, mt == MemoryType::PRESERVE ? mat.A.GetMemoryType() : mt);
195 A.CopyFrom(mat.A, nnz);
196
197 Rows = NULL;
198 #ifdef MFEM_USE_MEMALLOC
199 NodesMem = NULL;
200 #endif
201 }
202 else
203 {
204 #ifdef MFEM_USE_MEMALLOC
205 NodesMem = new RowNodeAlloc;
206 #endif
207 Rows = new RowNode *[height];
208 for (int i = 0; i < height; i++)
209 {
210 RowNode **node_pp = &Rows[i];
211 for (RowNode *node_p = mat.Rows[i]; node_p; node_p = node_p->Prev)
212 {
213 #ifdef MFEM_USE_MEMALLOC
214 RowNode *new_node_p = NodesMem->Alloc();
215 #else
216 RowNode *new_node_p = new RowNode;
217 #endif
218 new_node_p->Value = node_p->Value;
219 new_node_p->Column = node_p->Column;
220 *node_pp = new_node_p;
221 node_pp = &new_node_p->Prev;
222 }
223 *node_pp = NULL;
224 }
225
226 // We probably do not need to set the ownership flags here.
227 I.Reset(); I.SetHostPtrOwner(true);
228 J.Reset(); J.SetHostPtrOwner(true);
229 A.Reset(); A.SetHostPtrOwner(true);
230 }
231
232 current_row = -1;
233 ColPtrJ = NULL;
234 ColPtrNode = NULL;
235 At = NULL;
236 isSorted = mat.isSorted;
237
238 InitCuSparse();
239 }
240
SparseMatrix(const Vector & v)241 SparseMatrix::SparseMatrix(const Vector &v)
242 : AbstractSparseMatrix(v.Size(), v.Size())
243 , Rows(NULL)
244 , ColPtrJ(NULL)
245 , ColPtrNode(NULL)
246 , At(NULL)
247 , isSorted(true)
248 {
249 #ifdef MFEM_USE_MEMALLOC
250 NodesMem = NULL;
251 #endif
252 I.New(height + 1);
253 J.New(height);
254 A.New(height);
255
256 for (int i = 0; i <= height; i++)
257 {
258 I[i] = i;
259 }
260
261 for (int r=0; r<height; r++)
262 {
263 J[r] = r;
264 A[r] = v[r];
265 }
266
267 InitCuSparse();
268 }
269
operator =(const SparseMatrix & rhs)270 SparseMatrix& SparseMatrix::operator=(const SparseMatrix &rhs)
271 {
272 Clear();
273
274 SparseMatrix copy(rhs);
275 Swap(copy);
276
277 return *this;
278 }
279
MakeRef(const SparseMatrix & master)280 void SparseMatrix::MakeRef(const SparseMatrix &master)
281 {
282 MFEM_ASSERT(master.Finalized(), "'master' must be finalized");
283 Clear();
284 height = master.Height();
285 width = master.Width();
286 I = master.I; I.ClearOwnerFlags();
287 J = master.J; J.ClearOwnerFlags();
288 A = master.A; A.ClearOwnerFlags();
289 isSorted = master.isSorted;
290 }
291
SetEmpty()292 void SparseMatrix::SetEmpty()
293 {
294 height = width = 0;
295 I.Reset();
296 J.Reset();
297 A.Reset();
298 Rows = NULL;
299 current_row = -1;
300 ColPtrJ = NULL;
301 ColPtrNode = NULL;
302 At = NULL;
303 #ifdef MFEM_USE_MEMALLOC
304 NodesMem = NULL;
305 #endif
306 isSorted = false;
307
308 ClearCuSparse();
309 }
310
RowSize(const int i) const311 int SparseMatrix::RowSize(const int i) const
312 {
313 int gi = i;
314 if (gi < 0)
315 {
316 gi = -1-gi;
317 }
318
319 if (I)
320 {
321 return I[gi+1]-I[gi];
322 }
323
324 int s = 0;
325 RowNode *row = Rows[gi];
326 for ( ; row != NULL; row = row->Prev)
327 if (row->Value != 0.0)
328 {
329 s++;
330 }
331 return s;
332 }
333
MaxRowSize() const334 int SparseMatrix::MaxRowSize() const
335 {
336 int out=0;
337 int rowSize=0;
338 if (I)
339 {
340 for (int i=0; i < height; ++i)
341 {
342 rowSize = I[i+1]-I[i];
343 out = (out > rowSize) ? out : rowSize;
344 }
345 }
346 else
347 {
348 for (int i=0; i < height; ++i)
349 {
350 rowSize = RowSize(i);
351 out = (out > rowSize) ? out : rowSize;
352 }
353 }
354
355 return out;
356 }
357
GetRowColumns(const int row)358 int *SparseMatrix::GetRowColumns(const int row)
359 {
360 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
361
362 return J + I[row];
363 }
364
GetRowColumns(const int row) const365 const int *SparseMatrix::GetRowColumns(const int row) const
366 {
367 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
368
369 return J + I[row];
370 }
371
GetRowEntries(const int row)372 double *SparseMatrix::GetRowEntries(const int row)
373 {
374 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
375
376 return A + I[row];
377 }
378
GetRowEntries(const int row) const379 const double *SparseMatrix::GetRowEntries(const int row) const
380 {
381 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
382
383 return A + I[row];
384 }
385
SetWidth(int newWidth)386 void SparseMatrix::SetWidth(int newWidth)
387 {
388 if (newWidth == width)
389 {
390 // Nothing to be done here
391 return;
392 }
393 else if (newWidth == -1)
394 {
395 // Compute the actual width
396 width = ActualWidth();
397 // No need to reset the ColPtr, since the new ColPtr will be shorter.
398 }
399 else if (newWidth > width)
400 {
401 // We need to reset ColPtr, since now we may have additional columns.
402 if (Rows != NULL)
403 {
404 delete [] ColPtrNode;
405 ColPtrNode = static_cast<RowNode **>(NULL);
406 }
407 else
408 {
409 delete [] ColPtrJ;
410 ColPtrJ = static_cast<int *>(NULL);
411 }
412 width = newWidth;
413 }
414 else
415 {
416 // Check that the new width is bigger or equal to the actual width.
417 MFEM_ASSERT(newWidth >= ActualWidth(),
418 "The new width needs to be bigger or equal to the actual width");
419 width = newWidth;
420 }
421 }
422
423
SortColumnIndices()424 void SparseMatrix::SortColumnIndices()
425 {
426 MFEM_VERIFY(Finalized(), "Matrix is not Finalized!");
427
428 if (isSorted)
429 {
430 return;
431 }
432
433 const int * Ip=HostReadI();
434 HostReadWriteJ();
435 HostReadWriteData();
436
437 Array<Pair<int,double> > row;
438 for (int j = 0, i = 0; i < height; i++)
439 {
440 int end = Ip[i+1];
441 row.SetSize(end - j);
442 for (int k = 0; k < row.Size(); k++)
443 {
444 row[k].one = J[j+k];
445 row[k].two = A[j+k];
446 }
447 row.Sort();
448 for (int k = 0; k < row.Size(); k++, j++)
449 {
450 J[j] = row[k].one;
451 A[j] = row[k].two;
452 }
453 }
454 isSorted = true;
455 }
456
MoveDiagonalFirst()457 void SparseMatrix::MoveDiagonalFirst()
458 {
459 MFEM_VERIFY(Finalized(), "Matrix is not Finalized!");
460
461 for (int row = 0, end = 0; row < height; row++)
462 {
463 int start = end, j;
464 end = I[row+1];
465 for (j = start; true; j++)
466 {
467 MFEM_VERIFY(j < end, "diagonal entry not found in row = " << row);
468 if (J[j] == row) { break; }
469 }
470 const double diag = A[j];
471 for ( ; j > start; j--)
472 {
473 J[j] = J[j-1];
474 A[j] = A[j-1];
475 }
476 J[start] = row;
477 A[start] = diag;
478 }
479 }
480
Elem(int i,int j)481 double &SparseMatrix::Elem(int i, int j)
482 {
483 return operator()(i,j);
484 }
485
Elem(int i,int j) const486 const double &SparseMatrix::Elem(int i, int j) const
487 {
488 return operator()(i,j);
489 }
490
operator ()(int i,int j)491 double &SparseMatrix::operator()(int i, int j)
492 {
493 MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
494 "Trying to access element outside of the matrix. "
495 << "height = " << height << ", "
496 << "width = " << width << ", "
497 << "i = " << i << ", "
498 << "j = " << j);
499
500 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
501
502 for (int k = I[i], end = I[i+1]; k < end; k++)
503 {
504 if (J[k] == j)
505 {
506 return A[k];
507 }
508 }
509
510 MFEM_ABORT("Did not find i = " << i << ", j = " << j << " in matrix.");
511 return A[0];
512 }
513
operator ()(int i,int j) const514 const double &SparseMatrix::operator()(int i, int j) const
515 {
516 static const double zero = 0.0;
517
518 MFEM_ASSERT(i < height && i >= 0 && j < width && j >= 0,
519 "Trying to access element outside of the matrix. "
520 << "height = " << height << ", "
521 << "width = " << width << ", "
522 << "i = " << i << ", "
523 << "j = " << j);
524
525 if (Finalized())
526 {
527 for (int k = I[i], end = I[i+1]; k < end; k++)
528 {
529 if (J[k] == j)
530 {
531 return A[k];
532 }
533 }
534 }
535 else
536 {
537 for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
538 {
539 if (node_p->Column == j)
540 {
541 return node_p->Value;
542 }
543 }
544 }
545
546 return zero;
547 }
548
GetDiag(Vector & d) const549 void SparseMatrix::GetDiag(Vector & d) const
550 {
551 MFEM_VERIFY(height == width, "Matrix must be square, not height = "
552 << height << ", width = " << width);
553 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
554
555 d.SetSize(height);
556
557 const auto I = this->ReadI();
558 const auto J = this->ReadJ();
559 const auto A = this->ReadData();
560 auto dd = d.Write();
561
562 MFEM_FORALL(i, height,
563 {
564 const int begin = I[i];
565 const int end = I[i+1];
566 int j;
567 for (j = begin; j < end; j++)
568 {
569 if (J[j] == i)
570 {
571 dd[i] = A[j];
572 break;
573 }
574 }
575 if (j == end)
576 {
577 dd[i] = 0.;
578 }
579 });
580 }
581
582 /// Produces a DenseMatrix from a SparseMatrix
ToDenseMatrix() const583 DenseMatrix *SparseMatrix::ToDenseMatrix() const
584 {
585 int num_rows = this->Height();
586 int num_cols = this->Width();
587
588 DenseMatrix * B = new DenseMatrix(num_rows, num_cols);
589
590 this->ToDenseMatrix(*B);
591
592 return B;
593 }
594
595 /// Produces a DenseMatrix from a SparseMatrix
ToDenseMatrix(DenseMatrix & B) const596 void SparseMatrix::ToDenseMatrix(DenseMatrix & B) const
597 {
598 B.SetSize(height, width);
599 B = 0.0;
600
601 for (int r=0; r<height; r++)
602 {
603 const int * col = this->GetRowColumns(r);
604 const double * val = this->GetRowEntries(r);
605
606 for (int cj=0; cj<this->RowSize(r); cj++)
607 {
608 B(r, col[cj]) = val[cj];
609 }
610 }
611 }
612
Mult(const Vector & x,Vector & y) const613 void SparseMatrix::Mult(const Vector &x, Vector &y) const
614 {
615 if (Finalized()) { y.UseDevice(true); }
616 y = 0.0;
617 AddMult(x, y);
618 }
619
AddMult(const Vector & x,Vector & y,const double a) const620 void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
621 {
622 MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
623 << ") must match matrix width (" << width << ")");
624 MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
625 << ") must match matrix height (" << height << ")");
626
627 if (!Finalized())
628 {
629 const double *xp = x.HostRead();
630 double *yp = y.HostReadWrite();
631
632 // The matrix is not finalized, but multiplication is still possible
633 for (int i = 0; i < height; i++)
634 {
635 RowNode *row = Rows[i];
636 double b = 0.0;
637 for ( ; row != NULL; row = row->Prev)
638 {
639 b += row->Value * xp[row->Column];
640 }
641 *yp += a * b;
642 yp++;
643 }
644 return;
645 }
646
647 #ifndef MFEM_USE_LEGACY_OPENMP
648 const int height = this->height;
649 const int nnz = J.Capacity();
650 auto d_I = Read(I, height+1);
651 auto d_J = Read(J, nnz);
652 auto d_A = Read(A, nnz);
653 auto d_x = x.Read();
654 auto d_y = y.ReadWrite();
655
656 // Skip if matrix has no non-zeros
657 if (nnz == 0) {return;}
658 if (Device::Allows(Backend::CUDA_MASK) && useCuSparse)
659 {
660 #ifdef MFEM_USE_CUDA
661 const double alpha = a;
662 const double beta = 1.0;
663
664 // Setup descriptors
665 if (!initBuffers)
666 {
667 #if CUDA_VERSION >= 10010
668 // Setup matrix descriptor
669 cusparseCreateCsr(&matA_descr,Height(), Width(), J.Capacity(),
670 const_cast<int *>(d_I),
671 const_cast<int *>(d_J), const_cast<double *>(d_A), CUSPARSE_INDEX_32I,
672 CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);
673
674 // Create handles for input/output vectors
675 cusparseCreateDnVec(&vecX_descr, x.Size(), const_cast<double *>(d_x),
676 CUDA_R_64F);
677 cusparseCreateDnVec(&vecY_descr, y.Size(), d_y, CUDA_R_64F);
678 #else
679 cusparseCreateMatDescr(&matA_descr);
680 cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
681 cusparseSetMatType(matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL);
682
683 #endif
684
685 initBuffers = true;
686 }
687 // Allocate kernel space. Buffer is shared between different sparsemats
688 size_t newBufferSize = 0;
689
690 #if CUDA_VERSION >= 11020
691 cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
692 matA_descr,
693 vecX_descr, &beta, vecY_descr, CUDA_R_64F,
694 CUSPARSE_SPMV_CSR_ALG1, &newBufferSize);
695 #elif CUDA_VERSION >= 10010
696 cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
697 matA_descr,
698 vecX_descr, &beta, vecY_descr, CUDA_R_64F,
699 CUSPARSE_CSRMV_ALG1, &newBufferSize);
700 #endif
701
702 // Check if we need to resize
703 if (newBufferSize > bufferSize)
704 {
705 bufferSize = newBufferSize;
706 if (dBuffer != nullptr) { CuMemFree(dBuffer); }
707 CuMemAlloc(&dBuffer, bufferSize);
708 }
709
710 #if CUDA_VERSION >= 11020
711 // Update input/output vectors
712 cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
713 cusparseDnVecSetValues(vecY_descr, d_y);
714
715 // Y = alpha A * X + beta * Y
716 cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
717 vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG1, dBuffer);
718 #elif CUDA_VERSION >= 10010
719 // Update input/output vectors
720 cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
721 cusparseDnVecSetValues(vecY_descr, d_y);
722
723 // Y = alpha A * X + beta * Y
724 cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
725 vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_CSRMV_ALG1, dBuffer);
726 #else
727 cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
728 Height(), Width(), J.Capacity(),
729 &alpha, matA_descr,
730 const_cast<double *>(d_A), const_cast<int *>(d_I), const_cast<int *>(d_J),
731 const_cast<double *>(d_x), &beta, d_y);
732 #endif
733 #endif
734 }
735 else
736 {
737 // Native version
738 MFEM_FORALL(i, height,
739 {
740 double d = 0.0;
741 const int end = d_I[i+1];
742 for (int j = d_I[i]; j < end; j++)
743 {
744 d += d_A[j] * d_x[d_J[j]];
745 }
746 d_y[i] += a * d;
747 });
748
749 }
750
751 #else
752 const double *Ap = A, *xp = x.GetData();
753 double *yp = y.GetData();
754 const int *Jp = J, *Ip = I;
755
756 #pragma omp parallel for
757 for (int i = 0; i < height; i++)
758 {
759 double d = 0.0;
760 const int end = Ip[i+1];
761 for (int j = Ip[i]; j < end; j++)
762 {
763 d += Ap[j] * xp[Jp[j]];
764 }
765 yp[i] += a * d;
766 }
767 #endif
768 }
769
MultTranspose(const Vector & x,Vector & y) const770 void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
771 {
772 if (Finalized()) { y.UseDevice(true); }
773 y = 0.0;
774 AddMultTranspose(x, y);
775 }
776
AddMultTranspose(const Vector & x,Vector & y,const double a) const777 void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
778 const double a) const
779 {
780 MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
781 << ") must match matrix height (" << height << ")");
782 MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
783 << ") must match matrix width (" << width << ")");
784
785 if (!Finalized())
786 {
787 double *yp = y.GetData();
788 // The matrix is not finalized, but multiplication is still possible
789 for (int i = 0; i < height; i++)
790 {
791 RowNode *row = Rows[i];
792 double b = a * x(i);
793 for ( ; row != NULL; row = row->Prev)
794 {
795 yp[row->Column] += row->Value * b;
796 }
797 }
798 return;
799 }
800
801 if (At)
802 {
803 At->AddMult(x, y, a);
804 }
805 else
806 {
807 MFEM_VERIFY(Device::IsDisabled(), "transpose action on device is not "
808 "enabled; see BuildTranspose() for details.");
809 for (int i = 0; i < height; i++)
810 {
811 const double xi = a * x[i];
812 const int end = I[i+1];
813 for (int j = I[i]; j < end; j++)
814 {
815 const int Jj = J[j];
816 y[Jj] += A[j] * xi;
817 }
818 }
819 }
820 }
821
BuildTranspose() const822 void SparseMatrix::BuildTranspose() const
823 {
824 if (At == NULL)
825 {
826 At = Transpose(*this);
827 }
828 }
829
ResetTranspose() const830 void SparseMatrix::ResetTranspose() const
831 {
832 delete At;
833 At = NULL;
834 }
835
PartMult(const Array<int> & rows,const Vector & x,Vector & y) const836 void SparseMatrix::PartMult(
837 const Array<int> &rows, const Vector &x, Vector &y) const
838 {
839 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
840
841 const int n = rows.Size();
842 const int nnz = J.Capacity();
843 auto d_rows = rows.Read();
844 auto d_I = Read(I, height+1);
845 auto d_J = Read(J, nnz);
846 auto d_A = Read(A, nnz);
847 auto d_x = x.Read();
848 auto d_y = y.Write();
849 MFEM_FORALL(i, n,
850 {
851 const int r = d_rows[i];
852 const int end = d_I[r + 1];
853 double a = 0.0;
854 for (int j = d_I[r]; j < end; j++)
855 {
856 a += d_A[j] * d_x[d_J[j]];
857 }
858 d_y[r] = a;
859 });
860 }
861
PartAddMult(const Array<int> & rows,const Vector & x,Vector & y,const double a) const862 void SparseMatrix::PartAddMult(
863 const Array<int> &rows, const Vector &x, Vector &y, const double a) const
864 {
865 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
866
867 for (int i = 0; i < rows.Size(); i++)
868 {
869 int r = rows[i];
870 int end = I[r + 1];
871 double val = 0.0;
872 for (int j = I[r]; j < end; j++)
873 {
874 val += A[j] * x(J[j]);
875 }
876 y(r) += a * val;
877 }
878 }
879
BooleanMult(const Array<int> & x,Array<int> & y) const880 void SparseMatrix::BooleanMult(const Array<int> &x, Array<int> &y) const
881 {
882 MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
883 MFEM_ASSERT(x.Size() == Width(), "Input vector size (" << x.Size()
884 << ") must match matrix width (" << Width() << ")");
885
886 y.SetSize(Height(), Device::GetDeviceMemoryType());
887
888 const int height = Height();
889 const int nnz = J.Capacity();
890 auto d_I = Read(I, height+1);
891 auto d_J = Read(J, nnz);
892 auto d_x = Read(x.GetMemory(), x.Size());
893 auto d_y = Write(y.GetMemory(), y.Size());
894 MFEM_FORALL(i, height,
895 {
896 bool d_yi = false;
897 const int end = d_I[i+1];
898 for (int j = d_I[i]; j < end; j++)
899 {
900 if (d_x[d_J[j]])
901 {
902 d_yi = true;
903 break;
904 }
905 }
906 d_y[i] = d_yi;
907 });
908 }
909
BooleanMultTranspose(const Array<int> & x,Array<int> & y) const910 void SparseMatrix::BooleanMultTranspose(const Array<int> &x,
911 Array<int> &y) const
912 {
913 MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
914 MFEM_ASSERT(x.Size() == Height(), "Input vector size (" << x.Size()
915 << ") must match matrix height (" << Height() << ")");
916
917 y.SetSize(Width());
918 y = 0;
919
920 for (int i = 0; i < Height(); i++)
921 {
922 if (x[i])
923 {
924 int end = I[i+1];
925 for (int j = I[i]; j < end; j++)
926 {
927 y[J[j]] = x[i];
928 }
929 }
930 }
931 }
932
AbsMult(const Vector & x,Vector & y) const933 void SparseMatrix::AbsMult(const Vector &x, Vector &y) const
934 {
935 MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
936 << ") must match matrix width (" << width << ")");
937 MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
938 << ") must match matrix height (" << height << ")");
939
940 if (Finalized()) { y.UseDevice(true); }
941 y = 0.0;
942
943 if (!Finalized())
944 {
945 const double *xp = x.HostRead();
946 double *yp = y.HostReadWrite();
947
948 // The matrix is not finalized, but multiplication is still possible
949 for (int i = 0; i < height; i++)
950 {
951 RowNode *row = Rows[i];
952 double b = 0.0;
953 for ( ; row != NULL; row = row->Prev)
954 {
955 b += std::abs(row->Value) * xp[row->Column];
956 }
957 *yp += b;
958 yp++;
959 }
960 return;
961 }
962
963 const int height = this->height;
964 const int nnz = J.Capacity();
965 auto d_I = Read(I, height+1);
966 auto d_J = Read(J, nnz);
967 auto d_A = Read(A, nnz);
968 auto d_x = x.Read();
969 auto d_y = y.ReadWrite();
970 MFEM_FORALL(i, height,
971 {
972 double d = 0.0;
973 const int end = d_I[i+1];
974 for (int j = d_I[i]; j < end; j++)
975 {
976 d += std::abs(d_A[j]) * d_x[d_J[j]];
977 }
978 d_y[i] += d;
979 });
980 }
981
AbsMultTranspose(const Vector & x,Vector & y) const982 void SparseMatrix::AbsMultTranspose(const Vector &x, Vector &y) const
983 {
984 MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
985 << ") must match matrix height (" << height << ")");
986 MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
987 << ") must match matrix width (" << width << ")");
988
989 y = 0.0;
990
991 if (!Finalized())
992 {
993 double *yp = y.GetData();
994 // The matrix is not finalized, but multiplication is still possible
995 for (int i = 0; i < height; i++)
996 {
997 RowNode *row = Rows[i];
998 double b = x(i);
999 for ( ; row != NULL; row = row->Prev)
1000 {
1001 yp[row->Column] += fabs(row->Value) * b;
1002 }
1003 }
1004 return;
1005 }
1006
1007 if (At)
1008 {
1009 At->AbsMult(x, y);
1010 }
1011 else
1012 {
1013 MFEM_VERIFY(Device::IsDisabled(), "transpose action on device is not "
1014 "enabled; see BuildTranspose() for details.");
1015 for (int i = 0; i < height; i++)
1016 {
1017 const double xi = x[i];
1018 const int end = I[i+1];
1019 for (int j = I[i]; j < end; j++)
1020 {
1021 const int Jj = J[j];
1022 y[Jj] += std::abs(A[j]) * xi;
1023 }
1024 }
1025 }
1026 }
1027
InnerProduct(const Vector & x,const Vector & y) const1028 double SparseMatrix::InnerProduct(const Vector &x, const Vector &y) const
1029 {
1030 MFEM_ASSERT(x.Size() == Width(), "x.Size() = " << x.Size()
1031 << " must be equal to Width() = " << Width());
1032 MFEM_ASSERT(y.Size() == Height(), "y.Size() = " << y.Size()
1033 << " must be equal to Height() = " << Height());
1034
1035 x.HostRead();
1036 y.HostRead();
1037 if (Finalized())
1038 {
1039 const int nnz = J.Capacity();
1040 HostRead(I, height+1);
1041 HostRead(J, nnz);
1042 HostRead(A, nnz);
1043 }
1044
1045 double prod = 0.0;
1046 for (int i = 0; i < height; i++)
1047 {
1048 double a = 0.0;
1049 if (A)
1050 {
1051 for (int j = I[i], end = I[i+1]; j < end; j++)
1052 {
1053 a += A[j] * x(J[j]);
1054 }
1055 }
1056 else
1057 {
1058 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
1059 {
1060 a += np->Value * x(np->Column);
1061 }
1062 }
1063 prod += a * y(i);
1064 }
1065
1066 return prod;
1067 }
1068
GetRowSums(Vector & x) const1069 void SparseMatrix::GetRowSums(Vector &x) const
1070 {
1071 for (int i = 0; i < height; i++)
1072 {
1073 double a = 0.0;
1074 if (A)
1075 {
1076 for (int j = I[i], end = I[i+1]; j < end; j++)
1077 {
1078 a += A[j];
1079 }
1080 }
1081 else
1082 {
1083 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
1084 {
1085 a += np->Value;
1086 }
1087 }
1088 x(i) = a;
1089 }
1090 }
1091
GetRowNorml1(int irow) const1092 double SparseMatrix::GetRowNorml1(int irow) const
1093 {
1094 MFEM_VERIFY(irow < height,
1095 "row " << irow << " not in matrix with height " << height);
1096
1097 double a = 0.0;
1098 if (A)
1099 {
1100 for (int j = I[irow], end = I[irow+1]; j < end; j++)
1101 {
1102 a += fabs(A[j]);
1103 }
1104 }
1105 else
1106 {
1107 for (RowNode *np = Rows[irow]; np != NULL; np = np->Prev)
1108 {
1109 a += fabs(np->Value);
1110 }
1111 }
1112
1113 return a;
1114 }
1115
Threshold(double tol,bool fix_empty_rows)1116 void SparseMatrix::Threshold(double tol, bool fix_empty_rows)
1117 {
1118 MFEM_ASSERT(Finalized(), "Matrix must be finalized.");
1119 double atol;
1120 atol = std::abs(tol);
1121
1122 fix_empty_rows = height == width ? fix_empty_rows : false;
1123
1124 double *newA;
1125 int *newI, *newJ;
1126 int i, j, nz;
1127
1128 newI = Memory<int>(height+1);
1129 newI[0] = 0;
1130 for (i = 0, nz = 0; i < height; i++)
1131 {
1132 bool found = false;
1133 for (j = I[i]; j < I[i+1]; j++)
1134 if (std::abs(A[j]) > atol)
1135 {
1136 found = true;
1137 nz++;
1138 }
1139 if (fix_empty_rows && !found) { nz++; }
1140 newI[i+1] = nz;
1141 }
1142
1143 newJ = Memory<int>(nz);
1144 newA = Memory<double>(nz);
1145 // Assume we're sorted until we find out otherwise
1146 isSorted = true;
1147 for (i = 0, nz = 0; i < height; i++)
1148 {
1149 bool found = false;
1150 int lastCol = -1;
1151 for (j = I[i]; j < I[i+1]; j++)
1152 if (std::abs(A[j]) > atol)
1153 {
1154 found = true;
1155 newJ[nz] = J[j];
1156 newA[nz] = A[j];
1157 if ( lastCol > newJ[nz] )
1158 {
1159 isSorted = false;
1160 }
1161 lastCol = newJ[nz];
1162 nz++;
1163 }
1164 if (fix_empty_rows && !found)
1165 {
1166 newJ[nz] = i;
1167 newA[nz] = 0.0;
1168 nz++;
1169 }
1170 }
1171 Destroy();
1172 I.Wrap(newI, height+1, true);
1173 J.Wrap(newJ, I[height], true);
1174 A.Wrap(newA, I[height], true);
1175 }
1176
Finalize(int skip_zeros,bool fix_empty_rows)1177 void SparseMatrix::Finalize(int skip_zeros, bool fix_empty_rows)
1178 {
1179 int i, j, nr, nz;
1180 RowNode *aux;
1181
1182 if (Finalized())
1183 {
1184 return;
1185 }
1186
1187 delete [] ColPtrNode;
1188 ColPtrNode = NULL;
1189
1190 I.New(height+1);
1191 I[0] = 0;
1192 for (i = 1; i <= height; i++)
1193 {
1194 nr = 0;
1195 for (aux = Rows[i-1]; aux != NULL; aux = aux->Prev)
1196 {
1197 if (skip_zeros && aux->Value == 0.0)
1198 {
1199 if (skip_zeros == 2) { continue; }
1200 if ((i-1) != aux->Column) { continue; }
1201
1202 bool found = false;
1203 double found_val;
1204 for (RowNode *other = Rows[aux->Column]; other != NULL; other = other->Prev)
1205 {
1206 if (other->Column == (i-1))
1207 {
1208 found = true;
1209 found_val = other->Value;
1210 break;
1211 }
1212 }
1213 if (found && found_val == 0.0) { continue; }
1214
1215 }
1216 nr++;
1217 }
1218 if (fix_empty_rows && !nr) { nr = 1; }
1219 I[i] = I[i-1] + nr;
1220 }
1221
1222 nz = I[height];
1223 J.New(nz);
1224 A.New(nz);
1225 // Assume we're sorted until we find out otherwise
1226 isSorted = true;
1227 for (j = i = 0; i < height; i++)
1228 {
1229 int lastCol = -1;
1230 nr = 0;
1231 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
1232 {
1233 if (skip_zeros && aux->Value == 0.0)
1234 {
1235 if (skip_zeros == 2) { continue; }
1236 if (i != aux->Column) { continue; }
1237
1238 bool found = false;
1239 double found_val;
1240 for (RowNode *other = Rows[aux->Column]; other != NULL; other = other->Prev)
1241 {
1242 if (other->Column == i)
1243 {
1244 found = true;
1245 found_val = other->Value;
1246 break;
1247 }
1248 }
1249 if (found && found_val == 0.0) { continue; }
1250 }
1251
1252 J[j] = aux->Column;
1253 A[j] = aux->Value;
1254
1255 if ( lastCol > J[j] )
1256 {
1257 isSorted = false;
1258 }
1259 lastCol = J[j];
1260
1261 j++;
1262 nr++;
1263 }
1264 if (fix_empty_rows && !nr)
1265 {
1266 J[j] = i;
1267 A[j] = 1.0;
1268 j++;
1269 }
1270 }
1271
1272 #ifdef MFEM_USE_MEMALLOC
1273 delete NodesMem;
1274 NodesMem = NULL;
1275 #else
1276 for (i = 0; i < height; i++)
1277 {
1278 RowNode *node_p = Rows[i];
1279 while (node_p != NULL)
1280 {
1281 aux = node_p;
1282 node_p = node_p->Prev;
1283 delete aux;
1284 }
1285 }
1286 #endif
1287
1288 delete [] Rows;
1289 Rows = NULL;
1290 }
1291
GetBlocks(Array2D<SparseMatrix * > & blocks) const1292 void SparseMatrix::GetBlocks(Array2D<SparseMatrix *> &blocks) const
1293 {
1294 int br = blocks.NumRows(), bc = blocks.NumCols();
1295 int nr = (height + br - 1)/br, nc = (width + bc - 1)/bc;
1296
1297 for (int j = 0; j < bc; j++)
1298 {
1299 for (int i = 0; i < br; i++)
1300 {
1301 int *bI = Memory<int>(nr + 1);
1302 for (int k = 0; k <= nr; k++)
1303 {
1304 bI[k] = 0;
1305 }
1306 blocks(i,j) = new SparseMatrix(bI, NULL, NULL, nr, nc);
1307 }
1308 }
1309
1310 for (int gr = 0; gr < height; gr++)
1311 {
1312 int bi = gr/nr, i = gr%nr + 1;
1313 if (Finalized())
1314 {
1315 for (int j = I[gr]; j < I[gr+1]; j++)
1316 {
1317 if (A[j] != 0.0)
1318 {
1319 blocks(bi, J[j]/nc)->I[i]++;
1320 }
1321 }
1322 }
1323 else
1324 {
1325 for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
1326 {
1327 if (n_p->Value != 0.0)
1328 {
1329 blocks(bi, n_p->Column/nc)->I[i]++;
1330 }
1331 }
1332 }
1333 }
1334
1335 for (int j = 0; j < bc; j++)
1336 {
1337 for (int i = 0; i < br; i++)
1338 {
1339 SparseMatrix &b = *blocks(i,j);
1340 int nnz = 0, rs;
1341 for (int k = 1; k <= nr; k++)
1342 {
1343 rs = b.I[k], b.I[k] = nnz, nnz += rs;
1344 }
1345 b.J.New(nnz);
1346 b.A.New(nnz);
1347 }
1348 }
1349
1350 for (int gr = 0; gr < height; gr++)
1351 {
1352 int bi = gr/nr, i = gr%nr + 1;
1353 if (Finalized())
1354 {
1355 for (int j = I[gr]; j < I[gr+1]; j++)
1356 {
1357 if (A[j] != 0.0)
1358 {
1359 SparseMatrix &b = *blocks(bi, J[j]/nc);
1360 b.J[b.I[i]] = J[j] % nc;
1361 b.A[b.I[i]] = A[j];
1362 b.I[i]++;
1363 }
1364 }
1365 }
1366 else
1367 {
1368 for (RowNode *n_p = Rows[gr]; n_p != NULL; n_p = n_p->Prev)
1369 {
1370 if (n_p->Value != 0.0)
1371 {
1372 SparseMatrix &b = *blocks(bi, n_p->Column/nc);
1373 b.J[b.I[i]] = n_p->Column % nc;
1374 b.A[b.I[i]] = n_p->Value;
1375 b.I[i]++;
1376 }
1377 }
1378 }
1379 }
1380 }
1381
IsSymmetric() const1382 double SparseMatrix::IsSymmetric() const
1383 {
1384 if (height != width)
1385 {
1386 return infinity();
1387 }
1388
1389 double symm = 0.0;
1390 if (Empty())
1391 {
1392 // return 0.0;
1393 }
1394 else if (Finalized())
1395 {
1396 for (int i = 1; i < height; i++)
1397 {
1398 for (int j = I[i]; j < I[i+1]; j++)
1399 {
1400 if (J[j] < i)
1401 {
1402 symm = std::max(symm, std::abs(A[j]-(*this)(J[j],i)));
1403 }
1404 }
1405 }
1406 }
1407 else
1408 {
1409 for (int i = 0; i < height; i++)
1410 {
1411 for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
1412 {
1413 int col = node_p->Column;
1414 if (col < i)
1415 {
1416 symm = std::max(symm, std::abs(node_p->Value-(*this)(col,i)));
1417 }
1418 }
1419 }
1420 }
1421 return symm;
1422 }
1423
Symmetrize()1424 void SparseMatrix::Symmetrize()
1425 {
1426 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
1427
1428 int i, j;
1429 for (i = 1; i < height; i++)
1430 {
1431 for (j = I[i]; j < I[i+1]; j++)
1432 {
1433 if (J[j] < i)
1434 {
1435 A[j] += (*this)(J[j],i);
1436 A[j] *= 0.5;
1437 (*this)(J[j],i) = A[j];
1438 }
1439 }
1440 }
1441 }
1442
NumNonZeroElems() const1443 int SparseMatrix::NumNonZeroElems() const
1444 {
1445 if (A != NULL) // matrix is finalized
1446 {
1447 return I[height];
1448 }
1449 else
1450 {
1451 int nnz = 0;
1452
1453 for (int i = 0; i < height; i++)
1454 {
1455 for (RowNode *node_p = Rows[i]; node_p != NULL; node_p = node_p->Prev)
1456 {
1457 nnz++;
1458 }
1459 }
1460
1461 return nnz;
1462 }
1463 }
1464
MaxNorm() const1465 double SparseMatrix::MaxNorm() const
1466 {
1467 double m = 0.0;
1468
1469 if (A)
1470 {
1471 int nnz = I[height];
1472 for (int j = 0; j < nnz; j++)
1473 {
1474 m = std::max(m, std::abs(A[j]));
1475 }
1476 }
1477 else
1478 {
1479 for (int i = 0; i < height; i++)
1480 {
1481 for (RowNode *n_p = Rows[i]; n_p != NULL; n_p = n_p->Prev)
1482 {
1483 m = std::max(m, std::abs(n_p->Value));
1484 }
1485 }
1486 }
1487 return m;
1488 }
1489
CountSmallElems(double tol) const1490 int SparseMatrix::CountSmallElems(double tol) const
1491 {
1492 int counter = 0;
1493
1494 if (A)
1495 {
1496 const int nz = I[height];
1497 const double *Ap = A;
1498
1499 for (int i = 0; i < nz; i++)
1500 {
1501 counter += (std::abs(Ap[i]) <= tol);
1502 }
1503 }
1504 else
1505 {
1506 for (int i = 0; i < height; i++)
1507 {
1508 for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1509 {
1510 counter += (std::abs(aux->Value) <= tol);
1511 }
1512 }
1513 }
1514
1515 return counter;
1516 }
1517
CheckFinite() const1518 int SparseMatrix::CheckFinite() const
1519 {
1520 if (Empty())
1521 {
1522 return 0;
1523 }
1524 else if (Finalized())
1525 {
1526 return mfem::CheckFinite(A, I[height]);
1527 }
1528 else
1529 {
1530 int counter = 0;
1531 for (int i = 0; i < height; i++)
1532 {
1533 for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1534 {
1535 counter += !IsFinite(aux->Value);
1536 }
1537 }
1538 return counter;
1539 }
1540 }
1541
Inverse() const1542 MatrixInverse *SparseMatrix::Inverse() const
1543 {
1544 return NULL;
1545 }
1546
EliminateRow(int row,const double sol,Vector & rhs)1547 void SparseMatrix::EliminateRow(int row, const double sol, Vector &rhs)
1548 {
1549 RowNode *aux;
1550
1551 MFEM_ASSERT(row < height && row >= 0,
1552 "Row " << row << " not in matrix of height " << height);
1553
1554 MFEM_VERIFY(!Finalized(), "Matrix must NOT be finalized.");
1555
1556 for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
1557 {
1558 rhs(aux->Column) -= sol * aux->Value;
1559 aux->Value = 0.0;
1560 }
1561 }
1562
EliminateRow(int row,DiagonalPolicy dpolicy)1563 void SparseMatrix::EliminateRow(int row, DiagonalPolicy dpolicy)
1564 {
1565 RowNode *aux;
1566
1567 MFEM_ASSERT(row < height && row >= 0,
1568 "Row " << row << " not in matrix of height " << height);
1569 MFEM_ASSERT(dpolicy != DIAG_KEEP, "Diagonal policy must not be DIAG_KEEP");
1570 MFEM_ASSERT(dpolicy != DIAG_ONE || height == width,
1571 "if dpolicy == DIAG_ONE, matrix must be square, not height = "
1572 << height << ", width = " << width);
1573
1574 if (Rows == NULL)
1575 {
1576 for (int i=I[row]; i < I[row+1]; ++i)
1577 {
1578 A[i]=0.0;
1579 }
1580 }
1581 else
1582 {
1583 for (aux = Rows[row]; aux != NULL; aux = aux->Prev)
1584 {
1585 aux->Value = 0.0;
1586 }
1587 }
1588
1589 if (dpolicy == DIAG_ONE)
1590 {
1591 SearchRow(row, row) = 1.;
1592 }
1593 }
1594
EliminateCol(int col,DiagonalPolicy dpolicy)1595 void SparseMatrix::EliminateCol(int col, DiagonalPolicy dpolicy)
1596 {
1597 MFEM_ASSERT(col < width && col >= 0,
1598 "Col " << col << " not in matrix of width " << width);
1599 MFEM_ASSERT(dpolicy != DIAG_KEEP, "Diagonal policy must not be DIAG_KEEP");
1600 MFEM_ASSERT(dpolicy != DIAG_ONE || height == width,
1601 "if dpolicy == DIAG_ONE, matrix must be square, not height = "
1602 << height << ", width = " << width);
1603
1604 if (Rows == NULL)
1605 {
1606 const int nnz = I[height];
1607 for (int jpos = 0; jpos != nnz; ++jpos)
1608 {
1609 if (J[jpos] == col)
1610 {
1611 A[jpos] = 0.0;
1612 }
1613 }
1614 }
1615 else
1616 {
1617 for (int i = 0; i < height; i++)
1618 {
1619 for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1620 {
1621 if (aux->Column == col)
1622 {
1623 aux->Value = 0.0;
1624 }
1625 }
1626 }
1627 }
1628
1629 if (dpolicy == DIAG_ONE)
1630 {
1631 SearchRow(col, col) = 1.0;
1632 }
1633 }
1634
EliminateCols(const Array<int> & cols,const Vector * x,Vector * b)1635 void SparseMatrix::EliminateCols(const Array<int> &cols, const Vector *x,
1636 Vector *b)
1637 {
1638 if (Rows == NULL)
1639 {
1640 for (int i = 0; i < height; i++)
1641 {
1642 for (int jpos = I[i]; jpos != I[i+1]; ++jpos)
1643 {
1644 if (cols[ J[jpos]])
1645 {
1646 if (x && b)
1647 {
1648 (*b)(i) -= A[jpos] * (*x)( J[jpos] );
1649 }
1650 A[jpos] = 0.0;
1651 }
1652 }
1653 }
1654 }
1655 else
1656 {
1657 for (int i = 0; i < height; i++)
1658 {
1659 for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
1660 {
1661 if (cols[aux -> Column])
1662 {
1663 if (x && b)
1664 {
1665 (*b)(i) -= aux -> Value * (*x)(aux -> Column);
1666 }
1667 aux->Value = 0.0;
1668 }
1669 }
1670 }
1671 }
1672 }
1673
EliminateCols(const Array<int> & col_marker,SparseMatrix & Ae)1674 void SparseMatrix::EliminateCols(const Array<int> &col_marker, SparseMatrix &Ae)
1675 {
1676 if (Rows)
1677 {
1678 RowNode *nd;
1679 for (int row = 0; row < height; row++)
1680 {
1681 for (nd = Rows[row]; nd != NULL; nd = nd->Prev)
1682 {
1683 if (col_marker[nd->Column])
1684 {
1685 Ae.Add(row, nd->Column, nd->Value);
1686 nd->Value = 0.0;
1687 }
1688 }
1689 }
1690 }
1691 else
1692 {
1693 for (int row = 0; row < height; row++)
1694 {
1695 for (int j = I[row]; j < I[row+1]; j++)
1696 {
1697 if (col_marker[J[j]])
1698 {
1699 Ae.Add(row, J[j], A[j]);
1700 A[j] = 0.0;
1701 }
1702 }
1703 }
1704 }
1705 }
1706
1707
EliminateRowCol(int rc,const double sol,Vector & rhs,DiagonalPolicy dpolicy)1708 void SparseMatrix::EliminateRowCol(int rc, const double sol, Vector &rhs,
1709 DiagonalPolicy dpolicy)
1710 {
1711 MFEM_ASSERT(rc < height && rc >= 0,
1712 "Row " << rc << " not in matrix of height " << height);
1713
1714 if (Rows == NULL)
1715 {
1716 for (int j = I[rc]; j < I[rc+1]; j++)
1717 {
1718 const int col = J[j];
1719 if (col == rc)
1720 {
1721 switch (dpolicy)
1722 {
1723 case DIAG_KEEP:
1724 rhs(rc) = A[j] * sol;
1725 break;
1726 case DIAG_ONE:
1727 A[j] = 1.0;
1728 rhs(rc) = sol;
1729 break;
1730 case DIAG_ZERO:
1731 A[j] = 0.;
1732 rhs(rc) = 0.;
1733 break;
1734 default:
1735 mfem_error("SparseMatrix::EliminateRowCol () #2");
1736 break;
1737 }
1738 }
1739 else
1740 {
1741 A[j] = 0.0;
1742 for (int k = I[col]; 1; k++)
1743 {
1744 if (k == I[col+1])
1745 {
1746 mfem_error("SparseMatrix::EliminateRowCol () #3");
1747 }
1748 else if (J[k] == rc)
1749 {
1750 rhs(col) -= sol * A[k];
1751 A[k] = 0.0;
1752 break;
1753 }
1754 }
1755 }
1756 }
1757 }
1758 else
1759 {
1760 for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
1761 {
1762 const int col = aux->Column;
1763 if (col == rc)
1764 {
1765 switch (dpolicy)
1766 {
1767 case DIAG_KEEP:
1768 rhs(rc) = aux->Value * sol;
1769 break;
1770 case DIAG_ONE:
1771 aux->Value = 1.0;
1772 rhs(rc) = sol;
1773 break;
1774 case DIAG_ZERO:
1775 aux->Value = 0.;
1776 rhs(rc) = 0.;
1777 break;
1778 default:
1779 mfem_error("SparseMatrix::EliminateRowCol () #4");
1780 break;
1781 }
1782 }
1783 else
1784 {
1785 aux->Value = 0.0;
1786 for (RowNode *node = Rows[col]; 1; node = node->Prev)
1787 {
1788 if (node == NULL)
1789 {
1790 mfem_error("SparseMatrix::EliminateRowCol () #5");
1791 }
1792 else if (node->Column == rc)
1793 {
1794 rhs(col) -= sol * node->Value;
1795 node->Value = 0.0;
1796 break;
1797 }
1798 }
1799 }
1800 }
1801 }
1802 }
1803
EliminateRowColMultipleRHS(int rc,const Vector & sol,DenseMatrix & rhs,DiagonalPolicy dpolicy)1804 void SparseMatrix::EliminateRowColMultipleRHS(int rc, const Vector &sol,
1805 DenseMatrix &rhs,
1806 DiagonalPolicy dpolicy)
1807 {
1808 MFEM_ASSERT(rc < height && rc >= 0,
1809 "Row " << rc << " not in matrix of height " << height);
1810 MFEM_ASSERT(sol.Size() == rhs.Width(), "solution size (" << sol.Size()
1811 << ") must match rhs width (" << rhs.Width() << ")");
1812
1813 const int num_rhs = rhs.Width();
1814 if (Rows == NULL)
1815 {
1816 for (int j = I[rc]; j < I[rc+1]; j++)
1817 {
1818 const int col = J[j];
1819 if (col == rc)
1820 {
1821 switch (dpolicy)
1822 {
1823 case DIAG_KEEP:
1824 for (int r = 0; r < num_rhs; r++)
1825 {
1826 rhs(rc,r) = A[j] * sol(r);
1827 }
1828 break;
1829 case DIAG_ONE:
1830 A[j] = 1.0;
1831 for (int r = 0; r < num_rhs; r++)
1832 {
1833 rhs(rc,r) = sol(r);
1834 }
1835 break;
1836 case DIAG_ZERO:
1837 A[j] = 0.;
1838 for (int r = 0; r < num_rhs; r++)
1839 {
1840 rhs(rc,r) = 0.;
1841 }
1842 break;
1843 default:
1844 mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #3");
1845 break;
1846 }
1847 }
1848 else
1849 {
1850 A[j] = 0.0;
1851 for (int k = I[col]; 1; k++)
1852 {
1853 if (k == I[col+1])
1854 {
1855 mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #4");
1856 }
1857 else if (J[k] == rc)
1858 {
1859 for (int r = 0; r < num_rhs; r++)
1860 {
1861 rhs(col,r) -= sol(r) * A[k];
1862 }
1863 A[k] = 0.0;
1864 break;
1865 }
1866 }
1867 }
1868 }
1869 }
1870 else
1871 {
1872 for (RowNode *aux = Rows[rc]; aux != NULL; aux = aux->Prev)
1873 {
1874 const int col = aux->Column;
1875 if (col == rc)
1876 {
1877 switch (dpolicy)
1878 {
1879 case DIAG_KEEP:
1880 for (int r = 0; r < num_rhs; r++)
1881 {
1882 rhs(rc,r) = aux->Value * sol(r);
1883 }
1884 break;
1885 case DIAG_ONE:
1886 aux->Value = 1.0;
1887 for (int r = 0; r < num_rhs; r++)
1888 {
1889 rhs(rc,r) = sol(r);
1890 }
1891 break;
1892 case DIAG_ZERO:
1893 aux->Value = 0.;
1894 for (int r = 0; r < num_rhs; r++)
1895 {
1896 rhs(rc,r) = 0.;
1897 }
1898 break;
1899 default:
1900 mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #5");
1901 break;
1902 }
1903 }
1904 else
1905 {
1906 aux->Value = 0.0;
1907 for (RowNode *node = Rows[col]; 1; node = node->Prev)
1908 {
1909 if (node == NULL)
1910 {
1911 mfem_error("SparseMatrix::EliminateRowColMultipleRHS() #6");
1912 }
1913 else if (node->Column == rc)
1914 {
1915 for (int r = 0; r < num_rhs; r++)
1916 {
1917 rhs(col,r) -= sol(r) * node->Value;
1918 }
1919 node->Value = 0.0;
1920 break;
1921 }
1922 }
1923 }
1924 }
1925 }
1926 }
1927
EliminateRowCol(int rc,DiagonalPolicy dpolicy)1928 void SparseMatrix::EliminateRowCol(int rc, DiagonalPolicy dpolicy)
1929 {
1930 MFEM_ASSERT(rc < height && rc >= 0,
1931 "Row " << rc << " not in matrix of height " << height);
1932
1933 if (Rows == NULL)
1934 {
1935 const auto &I = this->I; // only use const access for I
1936 const auto &J = this->J; // only use const access for J
1937 for (int j = I[rc]; j < I[rc+1]; j++)
1938 {
1939 const int col = J[j];
1940 if (col == rc)
1941 {
1942 if (dpolicy == DIAG_ONE)
1943 {
1944 A[j] = 1.0;
1945 }
1946 else if (dpolicy == DIAG_ZERO)
1947 {
1948 A[j] = 0.0;
1949 }
1950 }
1951 else
1952 {
1953 A[j] = 0.0;
1954 for (int k = I[col]; 1; k++)
1955 {
1956 if (k == I[col+1])
1957 {
1958 mfem_error("SparseMatrix::EliminateRowCol() #2");
1959 }
1960 else if (J[k] == rc)
1961 {
1962 A[k] = 0.0;
1963 break;
1964 }
1965 }
1966 }
1967 }
1968 }
1969 else
1970 {
1971 RowNode *aux, *node;
1972
1973 for (aux = Rows[rc]; aux != NULL; aux = aux->Prev)
1974 {
1975 const int col = aux->Column;
1976 if (col == rc)
1977 {
1978 if (dpolicy == DIAG_ONE)
1979 {
1980 aux->Value = 1.0;
1981 }
1982 else if (dpolicy == DIAG_ZERO)
1983 {
1984 aux->Value = 0.;
1985 }
1986 }
1987 else
1988 {
1989 aux->Value = 0.0;
1990 for (node = Rows[col]; 1; node = node->Prev)
1991 {
1992 if (node == NULL)
1993 {
1994 mfem_error("SparseMatrix::EliminateRowCol() #3");
1995 }
1996 else if (node->Column == rc)
1997 {
1998 node->Value = 0.0;
1999 break;
2000 }
2001 }
2002 }
2003 }
2004 }
2005 }
2006
2007 // This is almost identical to EliminateRowCol(int, int), except for
2008 // the A[j] = value; and aux->Value = value; lines.
EliminateRowColDiag(int rc,double value)2009 void SparseMatrix::EliminateRowColDiag(int rc, double value)
2010 {
2011 MFEM_ASSERT(rc < height && rc >= 0,
2012 "Row " << rc << " not in matrix of height " << height);
2013
2014 if (Rows == NULL)
2015 {
2016 for (int j = I[rc]; j < I[rc+1]; j++)
2017 {
2018 const int col = J[j];
2019 if (col == rc)
2020 {
2021 A[j] = value;
2022 }
2023 else
2024 {
2025 A[j] = 0.0;
2026 for (int k = I[col]; 1; k++)
2027 {
2028 if (k == I[col+1])
2029 {
2030 mfem_error("SparseMatrix::EliminateRowCol() #2");
2031 }
2032 else if (J[k] == rc)
2033 {
2034 A[k] = 0.0;
2035 break;
2036 }
2037 }
2038 }
2039 }
2040 }
2041 else
2042 {
2043 RowNode *aux, *node;
2044
2045 for (aux = Rows[rc]; aux != NULL; aux = aux->Prev)
2046 {
2047 const int col = aux->Column;
2048 if (col == rc)
2049 {
2050 aux->Value = value;
2051 }
2052 else
2053 {
2054 aux->Value = 0.0;
2055 for (node = Rows[col]; 1; node = node->Prev)
2056 {
2057 if (node == NULL)
2058 {
2059 mfem_error("SparseMatrix::EliminateRowCol() #3");
2060 }
2061 else if (node->Column == rc)
2062 {
2063 node->Value = 0.0;
2064 break;
2065 }
2066 }
2067 }
2068 }
2069 }
2070 }
2071
EliminateRowCol(int rc,SparseMatrix & Ae,DiagonalPolicy dpolicy)2072 void SparseMatrix::EliminateRowCol(int rc, SparseMatrix &Ae,
2073 DiagonalPolicy dpolicy)
2074 {
2075 if (Rows)
2076 {
2077 RowNode *nd, *nd2;
2078 for (nd = Rows[rc]; nd != NULL; nd = nd->Prev)
2079 {
2080 const int col = nd->Column;
2081 if (col == rc)
2082 {
2083 switch (dpolicy)
2084 {
2085 case DIAG_ONE:
2086 Ae.Add(rc, rc, nd->Value - 1.0);
2087 nd->Value = 1.0;
2088 break;
2089 case DIAG_ZERO:
2090 Ae.Add(rc, rc, nd->Value);
2091 nd->Value = 0.;
2092 break;
2093 case DIAG_KEEP:
2094 break;
2095 default:
2096 mfem_error("SparseMatrix::EliminateRowCol #1");
2097 break;
2098 }
2099 }
2100 else
2101 {
2102 Ae.Add(rc, col, nd->Value);
2103 nd->Value = 0.0;
2104 for (nd2 = Rows[col]; 1; nd2 = nd2->Prev)
2105 {
2106 if (nd2 == NULL)
2107 {
2108 mfem_error("SparseMatrix::EliminateRowCol #2");
2109 }
2110 else if (nd2->Column == rc)
2111 {
2112 Ae.Add(col, rc, nd2->Value);
2113 nd2->Value = 0.0;
2114 break;
2115 }
2116 }
2117 }
2118 }
2119 }
2120 else
2121 {
2122 for (int j = I[rc]; j < I[rc+1]; j++)
2123 {
2124 const int col = J[j];
2125 if (col == rc)
2126 {
2127 switch (dpolicy)
2128 {
2129 case DIAG_ONE:
2130 Ae.Add(rc, rc, A[j] - 1.0);
2131 A[j] = 1.0;
2132 break;
2133 case DIAG_ZERO:
2134 Ae.Add(rc, rc, A[j]);
2135 A[j] = 0.;
2136 break;
2137 case DIAG_KEEP:
2138 break;
2139 default:
2140 mfem_error("SparseMatrix::EliminateRowCol #3");
2141 break;
2142 }
2143 }
2144 else
2145 {
2146 Ae.Add(rc, col, A[j]);
2147 A[j] = 0.0;
2148 for (int k = I[col]; true; k++)
2149 {
2150 if (k == I[col+1])
2151 {
2152 mfem_error("SparseMatrix::EliminateRowCol #4");
2153 }
2154 else if (J[k] == rc)
2155 {
2156 Ae.Add(col, rc, A[k]);
2157 A[k] = 0.0;
2158 break;
2159 }
2160 }
2161 }
2162 }
2163 }
2164 }
2165
SetDiagIdentity()2166 void SparseMatrix::SetDiagIdentity()
2167 {
2168 for (int i = 0; i < height; i++)
2169 {
2170 if (I[i+1] == I[i]+1 && fabs(A[I[i]]) < 1e-16)
2171 {
2172 A[I[i]] = 1.0;
2173 }
2174 }
2175 }
2176
EliminateZeroRows(const double threshold)2177 void SparseMatrix::EliminateZeroRows(const double threshold)
2178 {
2179 for (int i = 0; i < height; i++)
2180 {
2181 double zero = 0.0;
2182 for (int j = I[i]; j < I[i+1]; j++)
2183 {
2184 zero += fabs(A[j]);
2185 }
2186 if (zero <= threshold)
2187 {
2188 for (int j = I[i]; j < I[i+1]; j++)
2189 {
2190 A[j] = (J[j] == i) ? 1.0 : 0.0;
2191 }
2192 }
2193 }
2194 }
2195
Gauss_Seidel_forw(const Vector & x,Vector & y) const2196 void SparseMatrix::Gauss_Seidel_forw(const Vector &x, Vector &y) const
2197 {
2198 if (!Finalized())
2199 {
2200 double *yp = y.GetData();
2201 const double *xp = x.GetData();
2202 RowNode *diag_p, *n_p, **R = Rows;
2203
2204 const int s = height;
2205 for (int i = 0; i < s; i++)
2206 {
2207 double sum = 0.0;
2208 diag_p = NULL;
2209 for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
2210 {
2211 const int c = n_p->Column;
2212 if (c == i)
2213 {
2214 diag_p = n_p;
2215 }
2216 else
2217 {
2218 sum += n_p->Value * yp[c];
2219 }
2220 }
2221
2222 if (diag_p != NULL && diag_p->Value != 0.0)
2223 {
2224 yp[i] = (xp[i] - sum) / diag_p->Value;
2225 }
2226 else if (xp[i] == sum)
2227 {
2228 yp[i] = sum;
2229 }
2230 else
2231 {
2232 mfem_error("SparseMatrix::Gauss_Seidel_forw()");
2233 }
2234 }
2235 }
2236 else
2237 {
2238 const int s = height;
2239 const int nnz = J.Capacity();
2240 const int *Ip = HostRead(I, s+1);
2241 const int *Jp = HostRead(J, nnz);
2242 const double *Ap = HostRead(A, nnz);
2243 double *yp = y.HostReadWrite();
2244 const double *xp = x.HostRead();
2245
2246 for (int i = 0, j = Ip[0]; i < s; i++)
2247 {
2248 const int end = Ip[i+1];
2249 double sum = 0.0;
2250 int d = -1;
2251 for ( ; j < end; j++)
2252 {
2253 const int c = Jp[j];
2254 if (c == i)
2255 {
2256 d = j;
2257 }
2258 else
2259 {
2260 sum += Ap[j] * yp[c];
2261 }
2262 }
2263
2264 if (d >= 0 && Ap[d] != 0.0)
2265 {
2266 yp[i] = (xp[i] - sum) / Ap[d];
2267 }
2268 else if (xp[i] == sum)
2269 {
2270 yp[i] = sum;
2271 }
2272 else
2273 {
2274 mfem_error("SparseMatrix::Gauss_Seidel_forw(...) #2");
2275 }
2276 }
2277 }
2278 }
2279
Gauss_Seidel_back(const Vector & x,Vector & y) const2280 void SparseMatrix::Gauss_Seidel_back(const Vector &x, Vector &y) const
2281 {
2282 if (!Finalized())
2283 {
2284 double *yp = y.GetData();
2285 const double *xp = x.GetData();
2286 RowNode *diag_p, *n_p, **R = Rows;
2287
2288 for (int i = height-1; i >= 0; i--)
2289 {
2290 double sum = 0.;
2291 diag_p = NULL;
2292 for (n_p = R[i]; n_p != NULL; n_p = n_p->Prev)
2293 {
2294 const int c = n_p->Column;
2295 if (c == i)
2296 {
2297 diag_p = n_p;
2298 }
2299 else
2300 {
2301 sum += n_p->Value * yp[c];
2302 }
2303 }
2304
2305 if (diag_p != NULL && diag_p->Value != 0.0)
2306 {
2307 yp[i] = (xp[i] - sum) / diag_p->Value;
2308 }
2309 else if (xp[i] == sum)
2310 {
2311 yp[i] = sum;
2312 }
2313 else
2314 {
2315 mfem_error("SparseMatrix::Gauss_Seidel_back()");
2316 }
2317 }
2318 }
2319 else
2320 {
2321 const int s = height;
2322 const int nnz = J.Capacity();
2323 const int *Ip = HostRead(I, s+1);
2324 const int *Jp = HostRead(J, nnz);
2325 const double *Ap = HostRead(A, nnz);
2326 double *yp = y.HostReadWrite();
2327 const double *xp = x.HostRead();
2328
2329 for (int i = s-1, j = Ip[s]-1; i >= 0; i--)
2330 {
2331 const int beg = Ip[i];
2332 double sum = 0.;
2333 int d = -1;
2334 for ( ; j >= beg; j--)
2335 {
2336 const int c = Jp[j];
2337 if (c == i)
2338 {
2339 d = j;
2340 }
2341 else
2342 {
2343 sum += Ap[j] * yp[c];
2344 }
2345 }
2346
2347 if (d >= 0 && Ap[d] != 0.0)
2348 {
2349 yp[i] = (xp[i] - sum) / Ap[d];
2350 }
2351 else if (xp[i] == sum)
2352 {
2353 yp[i] = sum;
2354 }
2355 else
2356 {
2357 mfem_error("SparseMatrix::Gauss_Seidel_back(...) #2");
2358 }
2359 }
2360 }
2361 }
2362
GetJacobiScaling() const2363 double SparseMatrix::GetJacobiScaling() const
2364 {
2365 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2366
2367 double sc = 1.0;
2368 for (int i = 0; i < height; i++)
2369 {
2370 int d = -1;
2371 double norm = 0.0;
2372 for (int j = I[i]; j < I[i+1]; j++)
2373 {
2374 if (J[j] == i)
2375 {
2376 d = j;
2377 }
2378 norm += fabs(A[j]);
2379 }
2380 if (d >= 0 && A[d] != 0.0)
2381 {
2382 double a = 1.8 * fabs(A[d]) / norm;
2383 if (a < sc)
2384 {
2385 sc = a;
2386 }
2387 }
2388 else
2389 {
2390 mfem_error("SparseMatrix::GetJacobiScaling() #2");
2391 }
2392 }
2393 return sc;
2394 }
2395
Jacobi(const Vector & b,const Vector & x0,Vector & x1,double sc) const2396 void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
2397 double sc) const
2398 {
2399 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2400
2401 for (int i = 0; i < height; i++)
2402 {
2403 int d = -1;
2404 double sum = b(i);
2405 for (int j = I[i]; j < I[i+1]; j++)
2406 {
2407 if (J[j] == i)
2408 {
2409 d = j;
2410 }
2411 else
2412 {
2413 sum -= A[j] * x0(J[j]);
2414 }
2415 }
2416 if (d >= 0 && A[d] != 0.0)
2417 {
2418 x1(i) = sc * (sum / A[d]) + (1.0 - sc) * x0(i);
2419 }
2420 else
2421 {
2422 mfem_error("SparseMatrix::Jacobi(...) #2");
2423 }
2424 }
2425 }
2426
DiagScale(const Vector & b,Vector & x,double sc) const2427 void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
2428 {
2429 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2430
2431 const int H = height;
2432 const int nnz = J.Capacity();
2433 const bool use_dev = b.UseDevice() || x.UseDevice();
2434
2435 const auto Ap = Read(A, nnz, use_dev);
2436 const auto Ip = Read(I, height+1, use_dev);
2437 const auto Jp = Read(J, nnz, use_dev);
2438
2439 const auto bp = b.Read(use_dev);
2440 auto xp = x.Write(use_dev);
2441
2442 MFEM_FORALL_SWITCH(use_dev, i, H,
2443 {
2444 const int end = Ip[i+1];
2445 for (int j = Ip[i]; true; j++)
2446 {
2447 if (j == end)
2448 {
2449 MFEM_ABORT_KERNEL("Diagonal not found in SparseMatrix::DiagScale");
2450 }
2451 if (Jp[j] == i)
2452 {
2453 if (!(std::abs(Ap[j]) > 0.0))
2454 {
2455 MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
2456 }
2457 xp[i] = sc * bp[i] / Ap[j];
2458 break;
2459 }
2460 }
2461 });
2462 }
2463
Jacobi2(const Vector & b,const Vector & x0,Vector & x1,double sc) const2464 void SparseMatrix::Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
2465 double sc) const
2466 {
2467 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2468
2469 for (int i = 0; i < height; i++)
2470 {
2471 double resi = b(i), norm = 0.0;
2472 for (int j = I[i]; j < I[i+1]; j++)
2473 {
2474 resi -= A[j] * x0(J[j]);
2475 norm += fabs(A[j]);
2476 }
2477 if (norm > 0.0)
2478 {
2479 x1(i) = x0(i) + sc * resi / norm;
2480 }
2481 else
2482 {
2483 MFEM_ABORT("L1 norm of row " << i << " is zero.");
2484 }
2485 }
2486 }
2487
Jacobi3(const Vector & b,const Vector & x0,Vector & x1,double sc) const2488 void SparseMatrix::Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
2489 double sc) const
2490 {
2491 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
2492
2493 for (int i = 0; i < height; i++)
2494 {
2495 double resi = b(i), sum = 0.0;
2496 for (int j = I[i]; j < I[i+1]; j++)
2497 {
2498 resi -= A[j] * x0(J[j]);
2499 sum += A[j];
2500 }
2501 if (sum > 0.0)
2502 {
2503 x1(i) = x0(i) + sc * resi / sum;
2504 }
2505 else
2506 {
2507 MFEM_ABORT("sum of row " << i << " is zero.");
2508 }
2509 }
2510 }
2511
AddSubMatrix(const Array<int> & rows,const Array<int> & cols,const DenseMatrix & subm,int skip_zeros)2512 void SparseMatrix::AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
2513 const DenseMatrix &subm, int skip_zeros)
2514 {
2515 int i, j, gi, gj, s, t;
2516 double a;
2517
2518 if (Finalized())
2519 {
2520 HostReadI();
2521 HostReadJ();
2522 HostReadWriteData();
2523 }
2524
2525 for (i = 0; i < rows.Size(); i++)
2526 {
2527 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2528 else { s = 1; }
2529 MFEM_ASSERT(gi < height,
2530 "Trying to insert a row " << gi << " outside the matrix height "
2531 << height);
2532 SetColPtr(gi);
2533 for (j = 0; j < cols.Size(); j++)
2534 {
2535 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2536 else { t = s; }
2537 MFEM_ASSERT(gj < width,
2538 "Trying to insert a column " << gj << " outside the matrix width "
2539 << width);
2540 a = subm(i, j);
2541 if (skip_zeros && a == 0.0)
2542 {
2543 // Skip assembly of zero elements if either:
2544 // (i) user specified to skip zeros regardless of symmetry, or
2545 // (ii) symmetry is not broken.
2546 if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2547 {
2548 continue;
2549 }
2550 }
2551 if (t < 0) { a = -a; }
2552 _Add_(gj, a);
2553 }
2554 ClearColPtr();
2555 }
2556 }
2557
Set(const int i,const int j,const double A)2558 void SparseMatrix::Set(const int i, const int j, const double A)
2559 {
2560 double a = A;
2561 int gi, gj, s, t;
2562
2563 if ((gi=i) < 0) { gi = -1-gi, s = -1; }
2564 else { s = 1; }
2565 MFEM_ASSERT(gi < height,
2566 "Trying to set a row " << gi << " outside the matrix height "
2567 << height);
2568 if ((gj=j) < 0) { gj = -1-gj, t = -s; }
2569 else { t = s; }
2570 MFEM_ASSERT(gj < width,
2571 "Trying to set a column " << gj << " outside the matrix width "
2572 << width);
2573 if (t < 0) { a = -a; }
2574 _Set_(gi, gj, a);
2575 }
2576
Add(const int i,const int j,const double A)2577 void SparseMatrix::Add(const int i, const int j, const double A)
2578 {
2579 int gi, gj, s, t;
2580 double a = A;
2581
2582 if ((gi=i) < 0) { gi = -1-gi, s = -1; }
2583 else { s = 1; }
2584 MFEM_ASSERT(gi < height,
2585 "Trying to insert a row " << gi << " outside the matrix height "
2586 << height);
2587 if ((gj=j) < 0) { gj = -1-gj, t = -s; }
2588 else { t = s; }
2589 MFEM_ASSERT(gj < width,
2590 "Trying to insert a column " << gj << " outside the matrix width "
2591 << width);
2592 if (t < 0) { a = -a; }
2593 _Add_(gi, gj, a);
2594 }
2595
SetSubMatrix(const Array<int> & rows,const Array<int> & cols,const DenseMatrix & subm,int skip_zeros)2596 void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
2597 const DenseMatrix &subm, int skip_zeros)
2598 {
2599 int i, j, gi, gj, s, t;
2600 double a;
2601
2602 for (i = 0; i < rows.Size(); i++)
2603 {
2604 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2605 else { s = 1; }
2606 MFEM_ASSERT(gi < height,
2607 "Trying to set a row " << gi << " outside the matrix height "
2608 << height);
2609 SetColPtr(gi);
2610 for (j = 0; j < cols.Size(); j++)
2611 {
2612 a = subm(i, j);
2613 if (skip_zeros && a == 0.0)
2614 {
2615 // Skip assembly of zero elements if either:
2616 // (i) user specified to skip zeros regardless of symmetry, or
2617 // (ii) symmetry is not broken.
2618 if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2619 {
2620 continue;
2621 }
2622 }
2623 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2624 else { t = s; }
2625 MFEM_ASSERT(gj < width,
2626 "Trying to set a column " << gj << " outside the matrix width "
2627 << width);
2628 if (t < 0) { a = -a; }
2629 _Set_(gj, a);
2630 }
2631 ClearColPtr();
2632 }
2633 }
2634
SetSubMatrixTranspose(const Array<int> & rows,const Array<int> & cols,const DenseMatrix & subm,int skip_zeros)2635 void SparseMatrix::SetSubMatrixTranspose(const Array<int> &rows,
2636 const Array<int> &cols,
2637 const DenseMatrix &subm,
2638 int skip_zeros)
2639 {
2640 int i, j, gi, gj, s, t;
2641 double a;
2642
2643 for (i = 0; i < rows.Size(); i++)
2644 {
2645 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2646 else { s = 1; }
2647 MFEM_ASSERT(gi < height,
2648 "Trying to set a row " << gi << " outside the matrix height "
2649 << height);
2650 SetColPtr(gi);
2651 for (j = 0; j < cols.Size(); j++)
2652 {
2653 a = subm(j, i);
2654 if (skip_zeros && a == 0.0)
2655 {
2656 // Skip assembly of zero elements if either:
2657 // (i) user specified to skip zeros regardless of symmetry, or
2658 // (ii) symmetry is not broken.
2659 if (skip_zeros == 2 || &rows != &cols || subm(j, i) == 0.0)
2660 {
2661 continue;
2662 }
2663 }
2664 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2665 else { t = s; }
2666 MFEM_ASSERT(gj < width,
2667 "Trying to set a column " << gj << " outside the matrix width "
2668 << width);
2669 if (t < 0) { a = -a; }
2670 _Set_(gj, a);
2671 }
2672 ClearColPtr();
2673 }
2674 }
2675
GetSubMatrix(const Array<int> & rows,const Array<int> & cols,DenseMatrix & subm) const2676 void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
2677 DenseMatrix &subm) const
2678 {
2679 int i, j, gi, gj, s, t;
2680 double a;
2681
2682 for (i = 0; i < rows.Size(); i++)
2683 {
2684 if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
2685 else { s = 1; }
2686 MFEM_ASSERT(gi < height,
2687 "Trying to read a row " << gi << " outside the matrix height "
2688 << height);
2689 SetColPtr(gi);
2690 for (j = 0; j < cols.Size(); j++)
2691 {
2692 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2693 else { t = s; }
2694 MFEM_ASSERT(gj < width,
2695 "Trying to read a column " << gj << " outside the matrix width "
2696 << width);
2697 a = _Get_(gj);
2698 subm(i, j) = (t < 0) ? (-a) : (a);
2699 }
2700 ClearColPtr();
2701 }
2702 }
2703
RowIsEmpty(const int row) const2704 bool SparseMatrix::RowIsEmpty(const int row) const
2705 {
2706 int gi;
2707
2708 if ((gi=row) < 0)
2709 {
2710 gi = -1-gi;
2711 }
2712 MFEM_ASSERT(gi < height,
2713 "Trying to query a row " << gi << " outside the matrix height "
2714 << height);
2715 if (Rows)
2716 {
2717 return (Rows[gi] == NULL);
2718 }
2719 else
2720 {
2721 return (I[gi] == I[gi+1]);
2722 }
2723 }
2724
GetRow(const int row,Array<int> & cols,Vector & srow) const2725 int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
2726 {
2727 RowNode *n;
2728 int j, gi;
2729
2730 if ((gi=row) < 0) { gi = -1-gi; }
2731 MFEM_ASSERT(gi < height,
2732 "Trying to read a row " << gi << " outside the matrix height "
2733 << height);
2734 if (Rows)
2735 {
2736 for (n = Rows[gi], j = 0; n; n = n->Prev)
2737 {
2738 j++;
2739 }
2740 cols.SetSize(j);
2741 srow.SetSize(j);
2742 for (n = Rows[gi], j = 0; n; n = n->Prev, j++)
2743 {
2744 cols[j] = n->Column;
2745 srow(j) = n->Value;
2746 }
2747 if (row < 0)
2748 {
2749 srow.Neg();
2750 }
2751
2752 return 0;
2753 }
2754 else
2755 {
2756 j = I[gi];
2757 cols.MakeRef(const_cast<int*>((const int*)J) + j, I[gi+1]-j);
2758 srow.NewDataAndSize(
2759 const_cast<double*>((const double*)A) + j, cols.Size());
2760 MFEM_ASSERT(row >= 0, "Row not valid: " << row << ", height: " << height);
2761 return 1;
2762 }
2763 }
2764
SetRow(const int row,const Array<int> & cols,const Vector & srow)2765 void SparseMatrix::SetRow(const int row, const Array<int> &cols,
2766 const Vector &srow)
2767 {
2768 int gi, gj, s, t;
2769 double a;
2770
2771 if ((gi=row) < 0) { gi = -1-gi, s = -1; }
2772 else { s = 1; }
2773 MFEM_ASSERT(gi < height,
2774 "Trying to set a row " << gi << " outside the matrix height "
2775 << height);
2776
2777 if (!Finalized())
2778 {
2779 SetColPtr(gi);
2780 for (int j = 0; j < cols.Size(); j++)
2781 {
2782 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2783 else { t = s; }
2784 MFEM_ASSERT(gj < width,
2785 "Trying to set a column " << gj << " outside the matrix"
2786 " width " << width);
2787 a = srow(j);
2788 if (t < 0) { a = -a; }
2789 _Set_(gj, a);
2790 }
2791 ClearColPtr();
2792 }
2793 else
2794 {
2795 MFEM_ASSERT(cols.Size() == RowSize(gi), "");
2796 MFEM_ASSERT(cols.Size() == srow.Size(), "");
2797
2798 for (int i = I[gi], j = 0; j < cols.Size(); j++, i++)
2799 {
2800 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2801 else { t = s; }
2802 MFEM_ASSERT(gj < width,
2803 "Trying to set a column " << gj << " outside the matrix"
2804 " width " << width);
2805
2806 J[i] = gj;
2807 A[i] = srow[j] * t;
2808 }
2809 }
2810 }
2811
AddRow(const int row,const Array<int> & cols,const Vector & srow)2812 void SparseMatrix::AddRow(const int row, const Array<int> &cols,
2813 const Vector &srow)
2814 {
2815 int j, gi, gj, s, t;
2816 double a;
2817
2818 MFEM_VERIFY(!Finalized(), "Matrix must NOT be finalized.");
2819
2820 if ((gi=row) < 0) { gi = -1-gi, s = -1; }
2821 else { s = 1; }
2822 MFEM_ASSERT(gi < height,
2823 "Trying to insert a row " << gi << " outside the matrix height "
2824 << height);
2825 SetColPtr(gi);
2826 for (j = 0; j < cols.Size(); j++)
2827 {
2828 if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
2829 else { t = s; }
2830 MFEM_ASSERT(gj < width,
2831 "Trying to insert a column " << gj << " outside the matrix width "
2832 << width);
2833 a = srow(j);
2834 if (a == 0.0)
2835 {
2836 continue;
2837 }
2838 if (t < 0) { a = -a; }
2839 _Add_(gj, a);
2840 }
2841 ClearColPtr();
2842 }
2843
ScaleRow(const int row,const double scale)2844 void SparseMatrix::ScaleRow(const int row, const double scale)
2845 {
2846 int i;
2847
2848 if ((i=row) < 0)
2849 {
2850 i = -1-i;
2851 }
2852 if (Rows != NULL)
2853 {
2854 RowNode *aux;
2855
2856 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
2857 {
2858 aux -> Value *= scale;
2859 }
2860 }
2861 else
2862 {
2863 int j, end = I[i+1];
2864
2865 for (j = I[i]; j < end; j++)
2866 {
2867 A[j] *= scale;
2868 }
2869 }
2870 }
2871
ScaleRows(const Vector & sl)2872 void SparseMatrix::ScaleRows(const Vector & sl)
2873 {
2874 double scale;
2875 if (Rows != NULL)
2876 {
2877 RowNode *aux;
2878 for (int i=0; i < height; ++i)
2879 {
2880 scale = sl(i);
2881 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
2882 {
2883 aux -> Value *= scale;
2884 }
2885 }
2886 }
2887 else
2888 {
2889 int j, end;
2890
2891 for (int i=0; i < height; ++i)
2892 {
2893 end = I[i+1];
2894 scale = sl(i);
2895 for (j = I[i]; j < end; j++)
2896 {
2897 A[j] *= scale;
2898 }
2899 }
2900 }
2901 }
2902
ScaleColumns(const Vector & sr)2903 void SparseMatrix::ScaleColumns(const Vector & sr)
2904 {
2905 if (Rows != NULL)
2906 {
2907 RowNode *aux;
2908 for (int i=0; i < height; ++i)
2909 {
2910 for (aux = Rows[i]; aux != NULL; aux = aux -> Prev)
2911 {
2912 aux -> Value *= sr(aux->Column);
2913 }
2914 }
2915 }
2916 else
2917 {
2918 int j, end;
2919
2920 for (int i=0; i < height; ++i)
2921 {
2922 end = I[i+1];
2923 for (j = I[i]; j < end; j++)
2924 {
2925 A[j] *= sr(J[j]);
2926 }
2927 }
2928 }
2929 }
2930
operator +=(const SparseMatrix & B)2931 SparseMatrix &SparseMatrix::operator+=(const SparseMatrix &B)
2932 {
2933 MFEM_ASSERT(height == B.height && width == B.width,
2934 "Mismatch of this matrix size and rhs. This height = "
2935 << height << ", width = " << width << ", B.height = "
2936 << B.height << ", B.width = " << width);
2937
2938 for (int i = 0; i < height; i++)
2939 {
2940 SetColPtr(i);
2941 if (B.Rows)
2942 {
2943 for (RowNode *aux = B.Rows[i]; aux != NULL; aux = aux->Prev)
2944 {
2945 _Add_(aux->Column, aux->Value);
2946 }
2947 }
2948 else
2949 {
2950 for (int j = B.I[i]; j < B.I[i+1]; j++)
2951 {
2952 _Add_(B.J[j], B.A[j]);
2953 }
2954 }
2955 ClearColPtr();
2956 }
2957
2958 return (*this);
2959 }
2960
Add(const double a,const SparseMatrix & B)2961 void SparseMatrix::Add(const double a, const SparseMatrix &B)
2962 {
2963 for (int i = 0; i < height; i++)
2964 {
2965 B.SetColPtr(i);
2966 if (Rows)
2967 {
2968 for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
2969 {
2970 np->Value += a * B._Get_(np->Column);
2971 }
2972 }
2973 else
2974 {
2975 for (int j = I[i]; j < I[i+1]; j++)
2976 {
2977 A[j] += a * B._Get_(J[j]);
2978 }
2979 }
2980 B.ClearColPtr();
2981 }
2982 }
2983
operator =(double a)2984 SparseMatrix &SparseMatrix::operator=(double a)
2985 {
2986 if (Rows == NULL)
2987 {
2988 const int nnz = J.Capacity();
2989 double *h_A = HostWrite(A, nnz);
2990 for (int i = 0; i < nnz; i++)
2991 {
2992 h_A[i] = a;
2993 }
2994 }
2995 else
2996 {
2997 for (int i = 0; i < height; i++)
2998 {
2999 for (RowNode *node_p = Rows[i]; node_p != NULL;
3000 node_p = node_p -> Prev)
3001 {
3002 node_p -> Value = a;
3003 }
3004 }
3005 }
3006
3007 return (*this);
3008 }
3009
operator *=(double a)3010 SparseMatrix &SparseMatrix::operator*=(double a)
3011 {
3012 if (Rows == NULL)
3013 {
3014 for (int i = 0, nnz = I[height]; i < nnz; i++)
3015 {
3016 A[i] *= a;
3017 }
3018 }
3019 else
3020 {
3021 for (int i = 0; i < height; i++)
3022 {
3023 for (RowNode *node_p = Rows[i]; node_p != NULL;
3024 node_p = node_p -> Prev)
3025 {
3026 node_p -> Value *= a;
3027 }
3028 }
3029 }
3030
3031 return (*this);
3032 }
3033
Print(std::ostream & out,int width_) const3034 void SparseMatrix::Print(std::ostream & out, int width_) const
3035 {
3036 int i, j;
3037
3038 if (A == NULL)
3039 {
3040 RowNode *nd;
3041 for (i = 0; i < height; i++)
3042 {
3043 out << "[row " << i << "]\n";
3044 for (nd = Rows[i], j = 0; nd != NULL; nd = nd->Prev, j++)
3045 {
3046 out << " (" << nd->Column << "," << nd->Value << ")";
3047 if ( !((j+1) % width_) )
3048 {
3049 out << '\n';
3050 }
3051 }
3052 if (j % width_)
3053 {
3054 out << '\n';
3055 }
3056 }
3057 return;
3058 }
3059
3060 // HostRead forces synchronization
3061 HostReadI();
3062 HostReadJ();
3063 HostReadData();
3064 for (i = 0; i < height; i++)
3065 {
3066 out << "[row " << i << "]\n";
3067 for (j = I[i]; j < I[i+1]; j++)
3068 {
3069 out << " (" << J[j] << "," << A[j] << ")";
3070 if ( !((j+1-I[i]) % width_) )
3071 {
3072 out << '\n';
3073 }
3074 }
3075 if ((j-I[i]) % width_)
3076 {
3077 out << '\n';
3078 }
3079 }
3080 }
3081
PrintMatlab(std::ostream & out) const3082 void SparseMatrix::PrintMatlab(std::ostream & out) const
3083 {
3084 out << "% size " << height << " " << width << "\n";
3085 out << "% Non Zeros " << NumNonZeroElems() << "\n";
3086 int i, j;
3087 ios::fmtflags old_fmt = out.flags();
3088 out.setf(ios::scientific);
3089 std::streamsize old_prec = out.precision(14);
3090
3091 for (i = 0; i < height; i++)
3092 {
3093 for (j = I[i]; j < I[i+1]; j++)
3094 {
3095 out << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
3096 }
3097 }
3098 out.precision(old_prec);
3099 out.flags(old_fmt);
3100 }
3101
PrintMM(std::ostream & out) const3102 void SparseMatrix::PrintMM(std::ostream & out) const
3103 {
3104 int i, j;
3105 ios::fmtflags old_fmt = out.flags();
3106 out.setf(ios::scientific);
3107 std::streamsize old_prec = out.precision(14);
3108
3109 out << "%%MatrixMarket matrix coordinate real general" << '\n'
3110 << "% Generated by MFEM" << '\n';
3111
3112 out << height << " " << width << " " << NumNonZeroElems() << '\n';
3113 for (i = 0; i < height; i++)
3114 {
3115 for (j = I[i]; j < I[i+1]; j++)
3116 {
3117 out << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
3118 }
3119 }
3120 out.precision(old_prec);
3121 out.flags(old_fmt);
3122 }
3123
PrintCSR(std::ostream & out) const3124 void SparseMatrix::PrintCSR(std::ostream & out) const
3125 {
3126 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
3127
3128 int i;
3129
3130 out << height << '\n'; // number of rows
3131
3132 for (i = 0; i <= height; i++)
3133 {
3134 out << I[i]+1 << '\n';
3135 }
3136
3137 for (i = 0; i < I[height]; i++)
3138 {
3139 out << J[i]+1 << '\n';
3140 }
3141
3142 for (i = 0; i < I[height]; i++)
3143 {
3144 out << A[i] << '\n';
3145 }
3146 }
3147
PrintCSR2(std::ostream & out) const3148 void SparseMatrix::PrintCSR2(std::ostream & out) const
3149 {
3150 MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
3151
3152 int i;
3153
3154 out << height << '\n'; // number of rows
3155 out << width << '\n'; // number of columns
3156
3157 for (i = 0; i <= height; i++)
3158 {
3159 out << I[i] << '\n';
3160 }
3161
3162 for (i = 0; i < I[height]; i++)
3163 {
3164 out << J[i] << '\n';
3165 }
3166
3167 for (i = 0; i < I[height]; i++)
3168 {
3169 out << A[i] << '\n';
3170 }
3171 }
3172
PrintInfo(std::ostream & out) const3173 void SparseMatrix::PrintInfo(std::ostream &out) const
3174 {
3175 const double MiB = 1024.*1024;
3176 int nnz = NumNonZeroElems();
3177 double pz = 100./nnz;
3178 int nz = CountSmallElems(0.0);
3179 double max_norm = MaxNorm();
3180 double symm = IsSymmetric();
3181 int nnf = CheckFinite();
3182 int ns12 = CountSmallElems(1e-12*max_norm);
3183 int ns15 = CountSmallElems(1e-15*max_norm);
3184 int ns18 = CountSmallElems(1e-18*max_norm);
3185
3186 out <<
3187 "SparseMatrix statistics:\n"
3188 " Format : " <<
3189 (Empty() ? "(empty)" : (Finalized() ? "CSR" : "LIL")) << "\n"
3190 " Dimensions : " << height << " x " << width << "\n"
3191 " Number of entries (total) : " << nnz << "\n"
3192 " Number of entries (per row) : " << 1.*nnz/Height() << "\n"
3193 " Number of stored zeros : " << nz*pz << "% (" << nz << ")\n"
3194 " Number of Inf/Nan entries : " << nnf*pz << "% ("<< nnf << ")\n"
3195 " Norm, max |a_ij| : " << max_norm << "\n"
3196 " Symmetry, max |a_ij-a_ji| : " << symm << "\n"
3197 " Number of small entries:\n"
3198 " |a_ij| <= 1e-12*Norm : " << ns12*pz << "% (" << ns12 << ")\n"
3199 " |a_ij| <= 1e-15*Norm : " << ns15*pz << "% (" << ns15 << ")\n"
3200 " |a_ij| <= 1e-18*Norm : " << ns18*pz << "% (" << ns18 << ")\n";
3201 if (Finalized())
3202 {
3203 out << " Memory used by CSR : " <<
3204 (sizeof(int)*(height+1+nnz)+sizeof(double)*nnz)/MiB << " MiB\n";
3205 }
3206 if (Rows != NULL)
3207 {
3208 size_t used_mem = sizeof(RowNode*)*height;
3209 #ifdef MFEM_USE_MEMALLOC
3210 used_mem += NodesMem->MemoryUsage();
3211 #else
3212 for (int i = 0; i < height; i++)
3213 {
3214 for (RowNode *aux = Rows[i]; aux != NULL; aux = aux->Prev)
3215 {
3216 used_mem += sizeof(RowNode);
3217 }
3218 }
3219 #endif
3220 out << " Memory used by LIL : " << used_mem/MiB << " MiB\n";
3221 }
3222 }
3223
Destroy()3224 void SparseMatrix::Destroy()
3225 {
3226 I.Delete();
3227 J.Delete();
3228 A.Delete();
3229
3230 if (Rows != NULL)
3231 {
3232 #if !defined(MFEM_USE_MEMALLOC)
3233 for (int i = 0; i < height; i++)
3234 {
3235 RowNode *aux, *node_p = Rows[i];
3236 while (node_p != NULL)
3237 {
3238 aux = node_p;
3239 node_p = node_p->Prev;
3240 delete aux;
3241 }
3242 }
3243 #endif
3244 delete [] Rows;
3245 }
3246
3247 delete [] ColPtrJ;
3248 delete [] ColPtrNode;
3249 #ifdef MFEM_USE_MEMALLOC
3250 delete NodesMem;
3251 #endif
3252 delete At;
3253
3254 ClearCuSparse();
3255 }
3256
ActualWidth() const3257 int SparseMatrix::ActualWidth() const
3258 {
3259 int awidth = 0;
3260 if (A)
3261 {
3262 const int *start_j = J;
3263 const int *end_j = J + I[height];
3264 for (const int *jptr = start_j; jptr != end_j; ++jptr)
3265 {
3266 awidth = std::max(awidth, *jptr + 1);
3267 }
3268 }
3269 else
3270 {
3271 RowNode *aux;
3272 for (int i = 0; i < height; i++)
3273 {
3274 for (aux = Rows[i]; aux != NULL; aux = aux->Prev)
3275 {
3276 awidth = std::max(awidth, aux->Column + 1);
3277 }
3278 }
3279 }
3280 return awidth;
3281 }
3282
SparseMatrixFunction(SparseMatrix & S,double (* f)(double))3283 void SparseMatrixFunction (SparseMatrix & S, double (*f)(double))
3284 {
3285 int n = S.NumNonZeroElems();
3286 double * s = S.GetData();
3287
3288 for (int i = 0; i < n; i++)
3289 {
3290 s[i] = f(s[i]);
3291 }
3292 }
3293
Transpose(const SparseMatrix & A)3294 SparseMatrix *Transpose (const SparseMatrix &A)
3295 {
3296 MFEM_VERIFY(
3297 A.Finalized(),
3298 "Finalize must be called before Transpose. Use TransposeRowMatrix instead");
3299
3300 int i, j, end;
3301 const int *A_i, *A_j;
3302 int m, n, nnz, *At_i, *At_j;
3303 const double *A_data;
3304 double *At_data;
3305
3306 m = A.Height(); // number of rows of A
3307 n = A.Width(); // number of columns of A
3308 nnz = A.NumNonZeroElems();
3309 A_i = A.GetI();
3310 A_j = A.GetJ();
3311 A_data = A.GetData();
3312
3313 At_i = Memory<int>(n+1);
3314 At_j = Memory<int>(nnz);
3315 At_data = Memory<double>(nnz);
3316
3317 for (i = 0; i <= n; i++)
3318 {
3319 At_i[i] = 0;
3320 }
3321 for (i = 0; i < nnz; i++)
3322 {
3323 At_i[A_j[i]+1]++;
3324 }
3325 for (i = 1; i < n; i++)
3326 {
3327 At_i[i+1] += At_i[i];
3328 }
3329
3330 for (i = j = 0; i < m; i++)
3331 {
3332 end = A_i[i+1];
3333 for ( ; j < end; j++)
3334 {
3335 At_j[At_i[A_j[j]]] = i;
3336 At_data[At_i[A_j[j]]] = A_data[j];
3337 At_i[A_j[j]]++;
3338 }
3339 }
3340
3341 for (i = n; i > 0; i--)
3342 {
3343 At_i[i] = At_i[i-1];
3344 }
3345 At_i[0] = 0;
3346
3347 return new SparseMatrix(At_i, At_j, At_data, n, m);
3348 }
3349
TransposeAbstractSparseMatrix(const AbstractSparseMatrix & A,int useActualWidth)3350 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
3351 int useActualWidth)
3352 {
3353 int i, j;
3354 int m, n, nnz, *At_i, *At_j;
3355 double *At_data;
3356 Array<int> Acols;
3357 Vector Avals;
3358
3359 m = A.Height(); // number of rows of A
3360 if (useActualWidth)
3361 {
3362 n = 0;
3363 int tmp;
3364 for (i = 0; i < m; i++)
3365 {
3366 A.GetRow(i, Acols, Avals);
3367 if (Acols.Size())
3368 {
3369 tmp = Acols.Max();
3370 if (tmp > n)
3371 {
3372 n = tmp;
3373 }
3374 }
3375 }
3376 ++n;
3377 }
3378 else
3379 {
3380 n = A.Width(); // number of columns of A
3381 }
3382 nnz = A.NumNonZeroElems();
3383
3384 At_i = Memory<int>(n+1);
3385 At_j = Memory<int>(nnz);
3386 At_data = Memory<double>(nnz);
3387
3388 for (i = 0; i <= n; i++)
3389 {
3390 At_i[i] = 0;
3391 }
3392
3393 for (i = 0; i < m; i++)
3394 {
3395 A.GetRow(i, Acols, Avals);
3396 for (j = 0; j<Acols.Size(); ++j)
3397 {
3398 At_i[Acols[j]+1]++;
3399 }
3400 }
3401 for (i = 1; i < n; i++)
3402 {
3403 At_i[i+1] += At_i[i];
3404 }
3405
3406 for (i = 0; i < m; i++)
3407 {
3408 A.GetRow(i, Acols, Avals);
3409 for (j = 0; j<Acols.Size(); ++j)
3410 {
3411 At_j[At_i[Acols[j]]] = i;
3412 At_data[At_i[Acols[j]]] = Avals[j];
3413 At_i[Acols[j]]++;
3414 }
3415 }
3416
3417 for (i = n; i > 0; i--)
3418 {
3419 At_i[i] = At_i[i-1];
3420 }
3421 At_i[0] = 0;
3422
3423 return new SparseMatrix(At_i, At_j, At_data, n, m);
3424 }
3425
3426
Mult(const SparseMatrix & A,const SparseMatrix & B,SparseMatrix * OAB)3427 SparseMatrix *Mult (const SparseMatrix &A, const SparseMatrix &B,
3428 SparseMatrix *OAB)
3429 {
3430 int nrowsA, ncolsA, nrowsB, ncolsB;
3431 const int *A_i, *A_j, *B_i, *B_j;
3432 int *C_i, *C_j, *B_marker;
3433 const double *A_data, *B_data;
3434 double *C_data;
3435 int ia, ib, ic, ja, jb, num_nonzeros;
3436 int row_start, counter;
3437 double a_entry, b_entry;
3438 SparseMatrix *C;
3439
3440 nrowsA = A.Height();
3441 ncolsA = A.Width();
3442 nrowsB = B.Height();
3443 ncolsB = B.Width();
3444
3445 MFEM_VERIFY(ncolsA == nrowsB,
3446 "number of columns of A (" << ncolsA
3447 << ") must equal number of rows of B (" << nrowsB << ")");
3448
3449 A_i = A.HostReadI();
3450 A_j = A.HostReadJ();
3451 A_data = A.HostReadData();
3452 B_i = B.HostReadI();
3453 B_j = B.HostReadJ();
3454 B_data = B.HostReadData();
3455
3456 B_marker = new int[ncolsB];
3457
3458 for (ib = 0; ib < ncolsB; ib++)
3459 {
3460 B_marker[ib] = -1;
3461 }
3462
3463 if (OAB == NULL)
3464 {
3465 C_i = Memory<int>(nrowsA+1);
3466
3467 C_i[0] = num_nonzeros = 0;
3468 for (ic = 0; ic < nrowsA; ic++)
3469 {
3470 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3471 {
3472 ja = A_j[ia];
3473 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
3474 {
3475 jb = B_j[ib];
3476 if (B_marker[jb] != ic)
3477 {
3478 B_marker[jb] = ic;
3479 num_nonzeros++;
3480 }
3481 }
3482 }
3483 C_i[ic+1] = num_nonzeros;
3484 }
3485
3486 C_j = Memory<int>(num_nonzeros);
3487 C_data = Memory<double>(num_nonzeros);
3488
3489 C = new SparseMatrix(C_i, C_j, C_data, nrowsA, ncolsB);
3490
3491 for (ib = 0; ib < ncolsB; ib++)
3492 {
3493 B_marker[ib] = -1;
3494 }
3495 }
3496 else
3497 {
3498 C = OAB;
3499
3500 MFEM_VERIFY(nrowsA == C->Height() && ncolsB == C->Width(),
3501 "Input matrix sizes do not match output sizes"
3502 << " nrowsA = " << nrowsA
3503 << ", C->Height() = " << C->Height()
3504 << " ncolsB = " << ncolsB
3505 << ", C->Width() = " << C->Width());
3506
3507 // C_i = C->HostReadI(); // not used
3508 C_j = C->HostWriteJ();
3509 C_data = C->HostWriteData();
3510 }
3511
3512 counter = 0;
3513 for (ic = 0; ic < nrowsA; ic++)
3514 {
3515 // row_start = C_i[ic];
3516 row_start = counter;
3517 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3518 {
3519 ja = A_j[ia];
3520 a_entry = A_data[ia];
3521 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
3522 {
3523 jb = B_j[ib];
3524 b_entry = B_data[ib];
3525 if (B_marker[jb] < row_start)
3526 {
3527 B_marker[jb] = counter;
3528 if (OAB == NULL)
3529 {
3530 C_j[counter] = jb;
3531 }
3532 C_data[counter] = a_entry*b_entry;
3533 counter++;
3534 }
3535 else
3536 {
3537 C_data[B_marker[jb]] += a_entry*b_entry;
3538 }
3539 }
3540 }
3541 }
3542
3543 MFEM_VERIFY(
3544 OAB == NULL || counter == OAB->NumNonZeroElems(),
3545 "With pre-allocated output matrix, number of non-zeros ("
3546 << OAB->NumNonZeroElems()
3547 << ") did not match number of entries changed from matrix-matrix multiply, "
3548 << counter);
3549
3550 delete [] B_marker;
3551
3552 return C;
3553 }
3554
TransposeMult(const SparseMatrix & A,const SparseMatrix & B)3555 SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
3556 {
3557 SparseMatrix *At = Transpose(A);
3558 SparseMatrix *AtB = Mult(*At, B);
3559 delete At;
3560 return AtB;
3561 }
3562
MultAbstractSparseMatrix(const AbstractSparseMatrix & A,const AbstractSparseMatrix & B)3563 SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
3564 const AbstractSparseMatrix &B)
3565 {
3566 int nrowsA, ncolsA, nrowsB, ncolsB;
3567 int *C_i, *C_j, *B_marker;
3568 double *C_data;
3569 int ia, ib, ic, ja, jb, num_nonzeros;
3570 int row_start, counter;
3571 double a_entry, b_entry;
3572 SparseMatrix *C;
3573
3574 nrowsA = A.Height();
3575 ncolsA = A.Width();
3576 nrowsB = B.Height();
3577 ncolsB = B.Width();
3578
3579 MFEM_VERIFY(ncolsA == nrowsB,
3580 "number of columns of A (" << ncolsA
3581 << ") must equal number of rows of B (" << nrowsB << ")");
3582
3583 B_marker = new int[ncolsB];
3584
3585 for (ib = 0; ib < ncolsB; ib++)
3586 {
3587 B_marker[ib] = -1;
3588 }
3589
3590 C_i = Memory<int>(nrowsA+1);
3591
3592 C_i[0] = num_nonzeros = 0;
3593
3594 Array<int> colsA, colsB;
3595 Vector dataA, dataB;
3596 for (ic = 0; ic < nrowsA; ic++)
3597 {
3598 A.GetRow(ic, colsA, dataA);
3599 for (ia = 0; ia < colsA.Size(); ia++)
3600 {
3601 ja = colsA[ia];
3602 B.GetRow(ja, colsB, dataB);
3603 for (ib = 0; ib < colsB.Size(); ib++)
3604 {
3605 jb = colsB[ib];
3606 if (B_marker[jb] != ic)
3607 {
3608 B_marker[jb] = ic;
3609 num_nonzeros++;
3610 }
3611 }
3612 }
3613 C_i[ic+1] = num_nonzeros;
3614 }
3615
3616 C_j = Memory<int>(num_nonzeros);
3617 C_data = Memory<double>(num_nonzeros);
3618
3619 C = new SparseMatrix(C_i, C_j, C_data, nrowsA, ncolsB);
3620
3621 for (ib = 0; ib < ncolsB; ib++)
3622 {
3623 B_marker[ib] = -1;
3624 }
3625
3626 counter = 0;
3627 for (ic = 0; ic < nrowsA; ic++)
3628 {
3629 row_start = counter;
3630 A.GetRow(ic, colsA, dataA);
3631 for (ia = 0; ia < colsA.Size(); ia++)
3632 {
3633 ja = colsA[ia];
3634 a_entry = dataA[ia];
3635 B.GetRow(ja, colsB, dataB);
3636 for (ib = 0; ib < colsB.Size(); ib++)
3637 {
3638 jb = colsB[ib];
3639 b_entry = dataB[ib];
3640 if (B_marker[jb] < row_start)
3641 {
3642 B_marker[jb] = counter;
3643 C_j[counter] = jb;
3644 C_data[counter] = a_entry*b_entry;
3645 counter++;
3646 }
3647 else
3648 {
3649 C_data[B_marker[jb]] += a_entry*b_entry;
3650 }
3651 }
3652 }
3653 }
3654
3655 delete [] B_marker;
3656
3657 return C;
3658 }
3659
Mult(const SparseMatrix & A,DenseMatrix & B)3660 DenseMatrix *Mult (const SparseMatrix &A, DenseMatrix &B)
3661 {
3662 DenseMatrix *C = new DenseMatrix(A.Height(), B.Width());
3663 Vector columnB, columnC;
3664 for (int j = 0; j < B.Width(); ++j)
3665 {
3666 B.GetColumnReference(j, columnB);
3667 C->GetColumnReference(j, columnC);
3668 A.Mult(columnB, columnC);
3669 }
3670 return C;
3671 }
3672
RAP(const SparseMatrix & A,DenseMatrix & P)3673 DenseMatrix *RAP (const SparseMatrix &A, DenseMatrix &P)
3674 {
3675 DenseMatrix R (P, 't'); // R = P^T
3676 DenseMatrix *AP = Mult (A, P);
3677 DenseMatrix *RAP_ = new DenseMatrix(R.Height(), AP->Width());
3678 Mult (R, *AP, *RAP_);
3679 delete AP;
3680 return RAP_;
3681 }
3682
RAP(DenseMatrix & A,const SparseMatrix & P)3683 DenseMatrix *RAP(DenseMatrix &A, const SparseMatrix &P)
3684 {
3685 SparseMatrix *R = Transpose(P);
3686 DenseMatrix *RA = Mult(*R, A);
3687 DenseMatrix AtP(*RA, 't');
3688 delete RA;
3689 DenseMatrix *RAtP = Mult(*R, AtP);
3690 delete R;
3691 DenseMatrix * RAP_ = new DenseMatrix(*RAtP, 't');
3692 delete RAtP;
3693 return RAP_;
3694 }
3695
RAP(const SparseMatrix & A,const SparseMatrix & R,SparseMatrix * ORAP)3696 SparseMatrix *RAP (const SparseMatrix &A, const SparseMatrix &R,
3697 SparseMatrix *ORAP)
3698 {
3699 SparseMatrix *P = Transpose (R);
3700 SparseMatrix *AP = Mult (A, *P);
3701 delete P;
3702 SparseMatrix *RAP_ = Mult (R, *AP, ORAP);
3703 delete AP;
3704 return RAP_;
3705 }
3706
RAP(const SparseMatrix & Rt,const SparseMatrix & A,const SparseMatrix & P)3707 SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
3708 const SparseMatrix &P)
3709 {
3710 SparseMatrix * R = Transpose(Rt);
3711 SparseMatrix * RA = Mult(*R,A);
3712 delete R;
3713 SparseMatrix * out = Mult(*RA, P);
3714 delete RA;
3715 return out;
3716 }
3717
Mult_AtDA(const SparseMatrix & A,const Vector & D,SparseMatrix * OAtDA)3718 SparseMatrix *Mult_AtDA (const SparseMatrix &A, const Vector &D,
3719 SparseMatrix *OAtDA)
3720 {
3721 int i, At_nnz, *At_j;
3722 double *At_data;
3723
3724 SparseMatrix *At = Transpose (A);
3725 At_nnz = At -> NumNonZeroElems();
3726 At_j = At -> GetJ();
3727 At_data = At -> GetData();
3728 for (i = 0; i < At_nnz; i++)
3729 {
3730 At_data[i] *= D(At_j[i]);
3731 }
3732 SparseMatrix *AtDA = Mult (*At, A, OAtDA);
3733 delete At;
3734 return AtDA;
3735 }
3736
Add(double a,const SparseMatrix & A,double b,const SparseMatrix & B)3737 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
3738 const SparseMatrix & B)
3739 {
3740 int nrows = A.Height();
3741 int ncols = A.Width();
3742
3743 int * C_i = Memory<int>(nrows+1);
3744 int * C_j;
3745 double * C_data;
3746
3747 const int *A_i = A.GetI();
3748 const int *A_j = A.GetJ();
3749 const double *A_data = A.GetData();
3750
3751 const int *B_i = B.GetI();
3752 const int *B_j = B.GetJ();
3753 const double *B_data = B.GetData();
3754
3755 int * marker = new int[ncols];
3756 std::fill(marker, marker+ncols, -1);
3757
3758 int num_nonzeros = 0, jcol;
3759 C_i[0] = 0;
3760 for (int ic = 0; ic < nrows; ic++)
3761 {
3762 for (int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3763 {
3764 jcol = A_j[ia];
3765 marker[jcol] = ic;
3766 num_nonzeros++;
3767 }
3768 for (int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
3769 {
3770 jcol = B_j[ib];
3771 if (marker[jcol] != ic)
3772 {
3773 marker[jcol] = ic;
3774 num_nonzeros++;
3775 }
3776 }
3777 C_i[ic+1] = num_nonzeros;
3778 }
3779
3780 C_j = Memory<int>(num_nonzeros);
3781 C_data = Memory<double>(num_nonzeros);
3782
3783 for (int ia = 0; ia < ncols; ia++)
3784 {
3785 marker[ia] = -1;
3786 }
3787
3788 int pos = 0;
3789 for (int ic = 0; ic < nrows; ic++)
3790 {
3791 for (int ia = A_i[ic]; ia < A_i[ic+1]; ia++)
3792 {
3793 jcol = A_j[ia];
3794 C_j[pos] = jcol;
3795 C_data[pos] = a*A_data[ia];
3796 marker[jcol] = pos;
3797 pos++;
3798 }
3799 for (int ib = B_i[ic]; ib < B_i[ic+1]; ib++)
3800 {
3801 jcol = B_j[ib];
3802 if (marker[jcol] < C_i[ic])
3803 {
3804 C_j[pos] = jcol;
3805 C_data[pos] = b*B_data[ib];
3806 marker[jcol] = pos;
3807 pos++;
3808 }
3809 else
3810 {
3811 C_data[marker[jcol]] += b*B_data[ib];
3812 }
3813 }
3814 }
3815
3816 delete[] marker;
3817 return new SparseMatrix(C_i, C_j, C_data, nrows, ncols);
3818 }
3819
Add(const SparseMatrix & A,const SparseMatrix & B)3820 SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B)
3821 {
3822 return Add(1.,A,1.,B);
3823 }
3824
Add(Array<SparseMatrix * > & Ai)3825 SparseMatrix * Add(Array<SparseMatrix *> & Ai)
3826 {
3827 MFEM_ASSERT(Ai.Size() > 0, "invalid size Ai.Size() = " << Ai.Size());
3828
3829 SparseMatrix * accumulate = Ai[0];
3830 SparseMatrix * result = accumulate;
3831
3832 for (int i=1; i < Ai.Size(); ++i)
3833 {
3834 result = Add(*accumulate, *Ai[i]);
3835 if (i != 1)
3836 {
3837 delete accumulate;
3838 }
3839
3840 accumulate = result;
3841 }
3842
3843 return result;
3844 }
3845
3846 /// B += alpha * A
Add(const SparseMatrix & A,double alpha,DenseMatrix & B)3847 void Add(const SparseMatrix &A,
3848 double alpha, DenseMatrix &B)
3849 {
3850 for (int r = 0; r < B.Height(); r++)
3851 {
3852 const int * colA = A.GetRowColumns(r);
3853 const double * valA = A.GetRowEntries(r);
3854 for (int i=0; i<A.RowSize(r); i++)
3855 {
3856 B(r, colA[i]) += alpha * valA[i];
3857 }
3858 }
3859 }
3860
3861 /// Produces a block matrix with blocks A_{ij}*B
OuterProduct(const DenseMatrix & A,const DenseMatrix & B)3862 DenseMatrix *OuterProduct(const DenseMatrix &A, const DenseMatrix &B)
3863 {
3864 int mA = A.Height(), nA = A.Width();
3865 int mB = B.Height(), nB = B.Width();
3866
3867 DenseMatrix *C = new DenseMatrix(mA * mB, nA * nB);
3868 *C = 0.0;
3869 for (int i=0; i<mA; i++)
3870 {
3871 for (int j=0; j<nA; j++)
3872 {
3873 C->AddMatrix(A(i,j), B, i * mB, j * nB);
3874 }
3875 }
3876 return C;
3877 }
3878
3879 /// Produces a block matrix with blocks A_{ij}*B
OuterProduct(const DenseMatrix & A,const SparseMatrix & B)3880 SparseMatrix *OuterProduct(const DenseMatrix &A, const SparseMatrix &B)
3881 {
3882 int mA = A.Height(), nA = A.Width();
3883 int mB = B.Height(), nB = B.Width();
3884
3885 SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
3886
3887 for (int i=0; i<mA; i++)
3888 {
3889 for (int j=0; j<nA; j++)
3890 {
3891 for (int r=0; r<mB; r++)
3892 {
3893 const int * colB = B.GetRowColumns(r);
3894 const double * valB = B.GetRowEntries(r);
3895
3896 for (int cj=0; cj<B.RowSize(r); cj++)
3897 {
3898 C->Set(i * mB + r, j * nB + colB[cj], A(i,j) * valB[cj]);
3899 }
3900 }
3901 }
3902 }
3903 C->Finalize();
3904
3905 return C;
3906 }
3907
3908 /// Produces a block matrix with blocks A_{ij}*B
OuterProduct(const SparseMatrix & A,const DenseMatrix & B)3909 SparseMatrix *OuterProduct(const SparseMatrix &A, const DenseMatrix &B)
3910 {
3911 int mA = A.Height(), nA = A.Width();
3912 int mB = B.Height(), nB = B.Width();
3913
3914 SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
3915
3916 for (int r=0; r<mA; r++)
3917 {
3918 const int * colA = A.GetRowColumns(r);
3919 const double * valA = A.GetRowEntries(r);
3920
3921 for (int aj=0; aj<A.RowSize(r); aj++)
3922 {
3923 for (int i=0; i<mB; i++)
3924 {
3925 for (int j=0; j<nB; j++)
3926 {
3927 C->Set(r * mB + i, colA[aj] * nB + j, valA[aj] * B(i, j));
3928 }
3929 }
3930 }
3931 }
3932 C->Finalize();
3933
3934 return C;
3935 }
3936
3937 /// Produces a block matrix with blocks A_{ij}*B
OuterProduct(const SparseMatrix & A,const SparseMatrix & B)3938 SparseMatrix *OuterProduct(const SparseMatrix &A, const SparseMatrix &B)
3939 {
3940 int mA = A.Height(), nA = A.Width();
3941 int mB = B.Height(), nB = B.Width();
3942
3943 SparseMatrix *C = new SparseMatrix(mA * mB, nA * nB);
3944
3945 for (int ar=0; ar<mA; ar++)
3946 {
3947 const int * colA = A.GetRowColumns(ar);
3948 const double * valA = A.GetRowEntries(ar);
3949
3950 for (int aj=0; aj<A.RowSize(ar); aj++)
3951 {
3952 for (int br=0; br<mB; br++)
3953 {
3954 const int * colB = B.GetRowColumns(br);
3955 const double * valB = B.GetRowEntries(br);
3956
3957 for (int bj=0; bj<B.RowSize(br); bj++)
3958 {
3959 C->Set(ar * mB + br, colA[aj] * nB + colB[bj],
3960 valA[aj] * valB[bj]);
3961 }
3962 }
3963 }
3964 }
3965 C->Finalize();
3966
3967 return C;
3968 }
3969
Swap(SparseMatrix & other)3970 void SparseMatrix::Swap(SparseMatrix &other)
3971 {
3972 mfem::Swap(width, other.width);
3973 mfem::Swap(height, other.height);
3974 mfem::Swap(I, other.I);
3975 mfem::Swap(J, other.J);
3976 mfem::Swap(A, other.A);
3977 mfem::Swap(Rows, other.Rows);
3978 mfem::Swap(current_row, other.current_row);
3979 mfem::Swap(ColPtrJ, other.ColPtrJ);
3980 mfem::Swap(ColPtrNode, other.ColPtrNode);
3981 mfem::Swap(At, other.At);
3982
3983 #ifdef MFEM_USE_MEMALLOC
3984 mfem::Swap(NodesMem, other.NodesMem);
3985 #endif
3986
3987 mfem::Swap(isSorted, other.isSorted);
3988 }
3989
3990 }
3991