complex_matrices.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#include "complex_matrices.h"
27
28namespace oomph
29{
30 //============================================================================
31 /// Complete LU solve (overwrites RHS with solution). This is the
32 /// generic version which should not need to be over-written.
33 //============================================================================
34 void ComplexMatrixBase::solve(Vector<std::complex<double>>& rhs)
35 {
36#ifdef PARANOID
37 // Check Matrix is square
38 if (nrow() != ncol())
39 {
40 throw OomphLibError(
41 "This matrix is not square, the matrix MUST be square!",
42 OOMPH_CURRENT_FUNCTION,
43 OOMPH_EXCEPTION_LOCATION);
44 }
45 // Check that size of rhs = nrow()
46 if (rhs.size() != nrow())
47 {
48 std::ostringstream error_message_stream;
49 error_message_stream << "The rhs vector is not the right size. It is "
50 << rhs.size() << ", it should be " << nrow()
51 << std::endl;
52
53 throw OomphLibError(error_message_stream.str(),
54 OOMPH_CURRENT_FUNCTION,
55 OOMPH_EXCEPTION_LOCATION);
56 }
57#endif
58
59 // Call the LU decomposition
61
62 // Call the back substitution
63 lubksub(rhs);
64 }
65
66
67 //============================================================================
68 /// Complete LU solve (Nothing gets overwritten!). This generic
69 /// version should never need to be overwritten
70 //============================================================================
71 void ComplexMatrixBase::solve(const Vector<std::complex<double>>& rhs,
72 Vector<std::complex<double>>& soln)
73 {
74 // Set the solution vector equal to the rhs
75 // N.B. This won't work if we change to another vector format
76 soln = rhs;
77
78 // Overwrite the solution vector (rhs is unchanged)
79 solve(soln);
80 }
81
82
83 //=======================================================================
84 /// Delete the storage that has been allocated for the LU factors, if
85 /// the matrix data is not itself being overwritten.
86 //======================================================================
88 {
89 // Clean up the LU factor storage, if it has been allocated
90 // and it's not the same as the matrix storage
91 if ((!Overwrite_matrix_storage) && (LU_factors != 0))
92 {
93 // Delete the pointer to the LU factors
94 delete[] LU_factors;
95 // Null out the vector
96 LU_factors = 0;
97 }
98 }
99
100
101 //=======================================================================
102 /// Destructor clean up the LU factors that have been allocated
103 //======================================================================
105 {
106 // Delete the storage for the index vector
107 delete Index;
108
109 // Delete the pointer to the LU factors
110 delete[] LU_factors;
111
112 // Null out the vector
113 LU_factors = 0;
114 }
115
116 //============================================================================
117 /// LU decompose a matrix, over-writing it and recording the pivots
118 /// numbers in the Index vector.
119 /// Returns the sign of the determinant.
120 //============================================================================
122 {
123#ifdef PARANOID
124 // Check Matrix is square
125 if (N != M)
126 {
127 throw OomphLibError(
128 "This matrix is not square, the matrix MUST be square!",
129 OOMPH_CURRENT_FUNCTION,
130 OOMPH_EXCEPTION_LOCATION);
131 }
132#endif
133
134 // Small constant
135 const double small_number = 1.0e-20;
136
137 // If the Index vector has not already been allocated,
138 // allocated storage for the index of permutations
139 if (Index == 0)
140 {
141 Index = new Vector<long>;
142 }
143
144 // Resize the index vector to the correct length
145 Index->resize(N);
146
147 // Vector scaling stores the implicit scaling of each row
148 Vector<double> scaling(N);
149
150 // Integer to store the sign that must multiply the determinant as
151 // a consequence of the row/column interchanges
152 int signature = 1;
153
154 // Loop over rows to get implicit scaling information
155 for (unsigned long i = 0; i < N; i++)
156 {
157 double big = 0.0;
158 for (unsigned long j = 0; j < M; j++)
159 {
160 double tmp = std::abs((*this)(i, j));
161 if (tmp > big) big = tmp;
162 }
163 if (big == 0.0)
164 {
165 throw OomphLibError(
166 "Singular Matrix", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
167 }
168 // Save the scaling
169 scaling[i] = 1.0 / big;
170 }
171
172 // Firsly, we shall delete any previous LU storage.
173 // If the user calls this function twice without changing the matrix
174 // then it is their own inefficiency, not ours (this time).
176
177 // Now if we're saving the LU factors separately (the default).
179 {
180 // Allocate storage for the LU factors
181 // Assign space for the n rows
182 LU_factors = new std::complex<double>[N * M];
183
184 // Now we know that memory has been allocated, copy over
185 // the matrix values
186 for (unsigned long i = 0; i < (N * M); i++)
187 {
189 }
190 }
191 // Otherwise the pointer to the LU_factors is the same as the
192 // matrix data
193 else
194 {
196 }
197
198 // Loop over columns
199 for (unsigned long j = 0; j < M; j++)
200 {
201 // Initialise imax
202 unsigned long imax = 0;
203
204 for (unsigned long i = 0; i < j; i++)
205 {
206 std::complex<double> sum = LU_factors[M * i + j];
207 for (unsigned long k = 0; k < i; k++)
208 {
209 sum -= LU_factors[M * i + k] * LU_factors[M * k + j];
210 }
211 LU_factors[M * i + j] = sum;
212 }
213
214 // Initialise search for largest pivot element
215 double largest_entry = 0.0;
216 for (unsigned long i = j; i < N; i++)
217 {
218 std::complex<double> sum = LU_factors[M * i + j];
219 for (unsigned long k = 0; k < j; k++)
220 {
221 sum -= LU_factors[M * i + k] * LU_factors[M * k + j];
222 }
223 LU_factors[M * i + j] = sum;
224 // Set temporary
225 double tmp = scaling[i] * std::abs(sum);
226 if (tmp >= largest_entry)
227 {
228 largest_entry = tmp;
229 imax = i;
230 }
231 }
232
233 // Test to see if we need to interchange rows
234 if (j != imax)
235 {
236 for (unsigned long k = 0; k < N; k++)
237 {
238 std::complex<double> tmp = LU_factors[M * imax + k];
239 LU_factors[M * imax + k] = LU_factors[M * j + k];
240 LU_factors[M * j + k] = tmp;
241 }
242 // Change the parity of signature
243 signature = -signature;
244 // Interchange scale factor
245 scaling[imax] = scaling[j];
246 }
247
248 // Set the index
249 (*Index)[j] = imax;
250 if (LU_factors[M * j + j] == 0.0)
251 {
252 LU_factors[M * j + j] = small_number;
253 }
254 // Divide by pivot element
255 if (j != N - 1)
256 {
257 std::complex<double> tmp = 1.0 / LU_factors[M * j + j];
258 for (unsigned long i = j + 1; i < N; i++)
259 {
260 LU_factors[M * i + j] *= tmp;
261 }
262 }
263
264 } // End of loop over columns
265
266
267 // Now multiply all the diagonal terms together to get the determinant
268 // Note that we need to use the mantissa, exponent formulation to
269 // avoid underflow errors. This is way more tedious in complex arithmetic.
270 std::complex<double> determinant_mantissa(1.0, 0.0);
271 std::complex<int> determinant_exponent(0, 0);
272 for (unsigned i = 0; i < N; i++)
273 {
274 // Get the next entry in mantissa exponent form
275 std::complex<double> entry;
276 int iexp_real = 0, iexp_imag = 0;
277 entry =
278 std::complex<double>(frexp(LU_factors[M * i + i].real(), &iexp_real),
279 frexp(LU_factors[M * i + i].imag(), &iexp_imag));
280
281 // Now calculate the four terms that contribute to the real
282 // and imaginary parts
283 // i.e. A + Bi * C + Di = AC - BD + i(AD + BC)
284 double temp_mantissa[4];
285 int temp_exp[4];
286
287 // Get the first contribution to the real part, AC
288 temp_mantissa[0] = determinant_mantissa.real() * entry.real();
289 temp_exp[0] = determinant_exponent.real() + iexp_real;
290 // Get the second contribution to the real part, DB
291 temp_mantissa[1] = determinant_mantissa.imag() * entry.imag();
292 temp_exp[1] = determinant_exponent.imag() + iexp_imag;
293
294 // Get the first contribution to the imaginary part, AD
295 temp_mantissa[2] = determinant_mantissa.real() * entry.imag();
296 temp_exp[2] = determinant_exponent.real() + iexp_imag;
297 // Get the second contribution to the imaginary part, BC
298 temp_mantissa[3] = determinant_mantissa.imag() * entry.real();
299 temp_exp[3] = determinant_exponent.imag() + iexp_real;
300
301 // Align the exponents in the two terms for each sum (real or imaginary)
302 // We always align up to the larger exponent
303 for (unsigned i = 0; i < 3; i += 2)
304 {
305 // If the first exponent is smaller, move it up
306 if (temp_exp[i] < temp_exp[i + 1])
307 {
308 // The smaller term
309 int lower = temp_exp[i];
310 // Loop over the difference in the exponents,
311 // scaling down the mantissa each time
312 for (int index = lower; index < temp_exp[i + 1]; ++index)
313 {
314 temp_mantissa[i] /= 2.0;
315 ++temp_exp[i];
316 }
317 }
318 // Otherwise the second exponent is smaller
319 else
320 {
321 // The smaller term is the other
322 int lower = temp_exp[i + 1];
323 // Loop over the difference in the exponents,
324 // scaling down the mantissa each time
325 for (int index = lower; index < temp_exp[i]; ++index)
326 {
327 temp_mantissa[i + 1] /= 2.0;
328 ++temp_exp[i + 1];
329 }
330 }
331 }
332
333 // Now combine the terms AC - BD
334 // and Combine the terms AD + BC
335 determinant_mantissa =
336 std::complex<double>(temp_mantissa[0] - temp_mantissa[1],
337 temp_mantissa[2] + temp_mantissa[3]);
338 // The exponents will be the same, so pick one
339 determinant_exponent = std::complex<int>(temp_exp[0], temp_exp[2]);
340
341 // Finally normalise the result
342 determinant_mantissa =
343 std::complex<double>(frexp(determinant_mantissa.real(), &iexp_real),
344 frexp(determinant_mantissa.imag(), &iexp_imag));
345
346 int temp_real = determinant_exponent.real() + iexp_real;
347 int temp_imag = determinant_exponent.imag() + iexp_imag;
348
349 determinant_exponent = std::complex<int>(temp_real, temp_imag);
350 }
351
352 // Integer to store the sign of the determinant
353 int sign = 0;
354 // Find the sign of the determinant (left or right half-plane)
355 if (std::abs(determinant_mantissa) > 0.0)
356 {
357 sign = 1;
358 }
359 if (std::abs(determinant_mantissa) < 0.0)
360 {
361 sign = -1;
362 }
363
364 // Multiply the sign by the signature
365 sign *= signature;
366
367 // Return the sign of the determinant
368 return sign;
369 }
370
371 //============================================================================
372 /// Back substitute an LU decomposed matrix.
373 //============================================================================
374 void DenseComplexMatrix::lubksub(Vector<std::complex<double>>& rhs)
375 {
376#ifdef PARANOID
377 // Check Matrix is square
378 if (N != M)
379 {
380 throw OomphLibError(
381 "This matrix is not square, the matrix MUST be square!",
382 OOMPH_CURRENT_FUNCTION,
383 OOMPH_EXCEPTION_LOCATION);
384 }
385 // Check that size of rhs=nrow() and index=nrow() are correct!
386 if (rhs.size() != N)
387 {
388 std::ostringstream error_message_stream;
389 error_message_stream << "The rhs vector is not the right size. It is "
390 << rhs.size() << ", it should be " << N << std::endl;
391
392 throw OomphLibError(error_message_stream.str(),
393 OOMPH_CURRENT_FUNCTION,
394 OOMPH_EXCEPTION_LOCATION);
395 }
396 if (Index == 0)
397 {
398 throw OomphLibError(
399 "Index vector has not been allocated. Have you called ludecompse()\n",
400 OOMPH_CURRENT_FUNCTION,
401 OOMPH_EXCEPTION_LOCATION);
402 }
403 if (Index->size() != N)
404 {
405 std::ostringstream error_message_stream;
406 error_message_stream << "The Index vector is not the right size. It is "
407 << Index->size() << ", it should be " << N
408 << std::endl;
409
410 throw OomphLibError(error_message_stream.str(),
411 OOMPH_CURRENT_FUNCTION,
412 OOMPH_EXCEPTION_LOCATION);
413 }
414#endif
415
416 unsigned long ii = 0;
417 for (unsigned i = 0; i < N; i++)
418 {
419 unsigned long ip = (*Index)[i];
420 std::complex<double> sum = rhs[ip];
421 rhs[ip] = rhs[i];
422
423 if (ii != 0)
424 {
425 for (unsigned long j = ii - 1; j < i; j++)
426 {
427 sum -= LU_factors[M * i + j] * rhs[j];
428 }
429 }
430 else if (sum != 0.0)
431 {
432 ii = i + 1;
433 }
434 rhs[i] = sum;
435 }
436
437 // Now do the back substitution
438 for (long i = long(N) - 1; i >= 0; i--)
439 {
440 std::complex<double> sum = rhs[i];
441 for (long j = i + 1; j < long(N); j++)
442 sum -= LU_factors[M * i + j] * rhs[j];
443 rhs[i] = sum / LU_factors[M * i + i];
444 }
445 }
446
447 //============================================================================
448 /// Find the residual of Ax=b, i.e. r=b-Ax
449 //============================================================================
450 void DenseComplexMatrix::residual(const Vector<std::complex<double>>& x,
451 const Vector<std::complex<double>>& rhs,
452 Vector<std::complex<double>>& residual)
453 {
454#ifdef PARANOID
455 // Check Matrix is square
456 if (N != M)
457 {
458 throw OomphLibError(
459 "This matrix is not square, the matrix MUST be square!",
460 OOMPH_CURRENT_FUNCTION,
461 OOMPH_EXCEPTION_LOCATION);
462 }
463 // Check that size of rhs = nrow()
464 if (rhs.size() != N)
465 {
466 std::ostringstream error_message_stream;
467 error_message_stream << "The rhs vector is not the right size. It is "
468 << rhs.size() << ", it should be " << N << std::endl;
469
470 throw OomphLibError(error_message_stream.str(),
471 OOMPH_CURRENT_FUNCTION,
472 OOMPH_EXCEPTION_LOCATION);
473 }
474 // Check that the size of x is correct
475 if (x.size() != N)
476 {
477 std::ostringstream error_message_stream;
478 error_message_stream << "The x vector is not the right size. It is "
479 << x.size() << ", it should be " << N << std::endl;
480
481 throw OomphLibError(error_message_stream.str(),
482 OOMPH_CURRENT_FUNCTION,
483 OOMPH_EXCEPTION_LOCATION);
484 }
485#endif
486 // If size of residual is wrong, correct it!
487 if (residual.size() != N)
488 {
489 residual.resize(N);
490 }
491
492 // Multiply the matrix by the vector x in residual vector
493 for (unsigned long i = 0; i < N; i++)
494 {
495 residual[i] = rhs[i];
496 for (unsigned long j = 0; j < M; j++)
497 {
498 residual[i] -= Matrixdata[M * i + j] * x[j];
499 }
500 }
501 }
502
503
504 //============================================================================
505 /// Multiply the matrix by the vector x: soln=Ax
506 //============================================================================
507 void DenseComplexMatrix::multiply(const Vector<std::complex<double>>& x,
508 Vector<std::complex<double>>& soln)
509 {
510#ifdef PARANOID
511 // Check to see if x.size() = ncol().
512 if (x.size() != M)
513 {
514 std::ostringstream error_message_stream;
515 error_message_stream << "The x vector is not the right size. It is "
516 << x.size() << ", it should be " << M << std::endl;
517
518 throw OomphLibError(error_message_stream.str(),
519 OOMPH_CURRENT_FUNCTION,
520 OOMPH_EXCEPTION_LOCATION);
521 }
522#endif
523
524 if (soln.size() != N)
525 {
526 // Resize and initialize the solution vector
527 soln.resize(N);
528 }
529
530 // Multiply the matrix A, by the vector x
531 for (unsigned long i = 0; i < N; i++)
532 {
533 soln[i] = 0.0;
534 for (unsigned long j = 0; j < M; j++)
535 {
536 soln[i] += Matrixdata[M * i + j] * x[j];
537 }
538 }
539 }
540
541
542 //=================================================================
543 /// Multiply the transposed matrix by the vector x: soln=A^T x
544 //=================================================================
546 const Vector<std::complex<double>>& x, Vector<std::complex<double>>& soln)
547 {
548#ifdef PARANOID
549 // Check to see x.size() = nrow()
550 if (x.size() != N)
551 {
552 std::ostringstream error_message_stream;
553 error_message_stream << "The x vector is not the right size. It is "
554 << x.size() << ", it should be " << N << std::endl;
555
556 throw OomphLibError(error_message_stream.str(),
557 OOMPH_CURRENT_FUNCTION,
558 OOMPH_EXCEPTION_LOCATION);
559 }
560#endif
561
562 if (soln.size() != M)
563 {
564 // Resize and initialize the solution vector
565 soln.resize(M);
566 }
567
568 // Initialise the solution
569 for (unsigned long i = 0; i < M; i++)
570 {
571 soln[i] = 0.0;
572 }
573
574
575 // Matrix vector product
576 for (unsigned long i = 0; i < N; i++)
577 {
578 for (unsigned long j = 0; j < M; j++)
579 {
580 soln[j] += Matrixdata[N * i + j] * x[i];
581 }
582 }
583 }
584
585 /// /////////////////////////////////////////////////////////////////////
586 /// /////////////////////////////////////////////////////////////////////
587 /// /////////////////////////////////////////////////////////////////////
588
589
590 //===================================================================
591 // Interface to SuperLU wrapper
592 //===================================================================
593 extern "C"
594 {
596 int*,
597 int*,
598 int*,
599 std::complex<double>*,
600 int*,
601 int*,
602 std::complex<double>*,
603 int*,
604 int*,
605 int*,
606 void*,
607 int*);
608 }
609
610
611 //===================================================================
612 /// Perform LU decomposition. Return the sign of the determinant
613 //===================================================================
615 {
616#ifdef PARANOID
617 if (N != M)
618 {
619 std::ostringstream error_message_stream;
620 error_message_stream << "Can only solve for quadratic matrices\n"
621 << "N, M " << N << " " << M << std::endl;
622
623 throw OomphLibError(error_message_stream.str(),
624 OOMPH_CURRENT_FUNCTION,
625 OOMPH_EXCEPTION_LOCATION);
626 }
627#endif
628
629 // SuperLU expects compressed column format: Set flag for
630 // transpose (0/1) = (false/true)
631 int transpose = 0;
632
633 // Doc (0/1) = (false/true)
634 int doc = 0;
636
637 // Cast to integers for stupid SuperLU
638 int n_aux = (int)N;
639 int nnz_aux = (int)Nnz;
640
641 // Integer to hold the sign of the determinant
642 int sign = 0;
643
644 // Just do the lu decompose phase (i=1)
645 int i = 1;
646 sign = superlu_complex(&i,
647 &n_aux,
648 &nnz_aux,
649 0,
650 Value,
651 Row_index,
653 0,
654 &n_aux,
655 &transpose,
656 &doc,
657 &F_factors,
658 &Info);
659
660 // Return the sign of the determinant
661 return sign;
662 }
663
664
665 //===================================================================
666 /// Do the backsubstitution
667 //===================================================================
668 void CCComplexMatrix::lubksub(Vector<std::complex<double>>& rhs)
669 {
670#ifdef PARANOID
671 if (N != rhs.size())
672 {
673 std::ostringstream error_message_stream;
674 error_message_stream << "Size of RHS doesn't match matrix size\n"
675 << "N, rhs.size() " << N << " " << rhs.size()
676 << std::endl;
677
678 throw OomphLibError(error_message_stream.str(),
679 OOMPH_CURRENT_FUNCTION,
680 OOMPH_EXCEPTION_LOCATION);
681 }
682 if (N != M)
683 {
684 std::ostringstream error_message_stream;
685 error_message_stream << "Can only solve for quadratic matrices\n"
686 << "N, M " << N << " " << M << std::endl;
687
688 throw OomphLibError(error_message_stream.str(),
689 OOMPH_CURRENT_FUNCTION,
690 OOMPH_EXCEPTION_LOCATION);
691 }
692#endif
693
694 /// RHS vector
695 std::complex<double>* b = new std::complex<double>[N];
696
697 // Copy across
698 for (unsigned long i = 0; i < N; i++)
699 {
700 b[i] = rhs[i];
701 }
702
703 // SuperLU expects compressed column format: Set flag for
704 // transpose (0/1) = (false/true)
705 int transpose = 0;
706
707 // Doc (0/1) = (false/true)
708 int doc = 0;
710
711 // Number of RHSs
712 int nrhs = 1;
713
714
715 // Cast to integers for stupid SuperLU
716 int n_aux = (int)N;
717 int nnz_aux = (int)Nnz;
718
719 // SuperLU in three steps (lu decompose, back subst, cleanup)
720 // Do the solve phase
721 int i = 2;
723 &n_aux,
724 &nnz_aux,
725 &nrhs,
726 Value,
727 Row_index,
729 b,
730 &n_aux,
731 &transpose,
732 &doc,
733 &F_factors,
734 &Info);
735
736 // Copy b into rhs vector
737 for (unsigned long i = 0; i < N; i++)
738 {
739 rhs[i] = b[i];
740 }
741
742 // Cleanup
743 delete[] b;
744 }
745
746
747 //===================================================================
748 /// Cleanup storage
749 //===================================================================
751 {
752 // Only bother if we've done an LU solve
753 if (F_factors != 0)
754 {
755 // SuperLU expects compressed column format: Set flag for
756 // transpose (0/1) = (false/true)
757 int transpose = 0;
758
759 // Doc (0/1) = (false/true)
760 int doc = 0;
762
763
764 // Cast to integers for stupid SuperLU
765 int n_aux = (int)N;
766 int nnz_aux = (int)Nnz;
767
768 // SuperLU in three steps (lu decompose, back subst, cleanup)
769 // Flag to indicate which solve step to do (1, 2 or 3)
770 int i = 3;
772 &n_aux,
773 &nnz_aux,
774 0,
775 Value,
776 Row_index,
778 0,
779 &n_aux,
780 &transpose,
781 &doc,
782 &F_factors,
783 &Info);
784 }
785 }
786
787
788 //===================================================================
789 /// Work out residual vector r = b-Ax for candidate solution x
790 //===================================================================
791 void CCComplexMatrix::residual(const Vector<std::complex<double>>& x,
792 const Vector<std::complex<double>>& rhs,
793 Vector<std::complex<double>>& residual)
794 {
795#ifdef PARANOID
796 // Check Matrix is square
797 if (N != M)
798 {
799 throw OomphLibError(
800 "This matrix is not square, the matrix MUST be square!",
801 OOMPH_CURRENT_FUNCTION,
802 OOMPH_EXCEPTION_LOCATION);
803 }
804 // Check that size of rhs = nrow()
805 if (rhs.size() != N)
806 {
807 std::ostringstream error_message_stream;
808 error_message_stream << "The rhs vector is not the right size. It is "
809 << rhs.size() << ", it should be " << N << std::endl;
810
811 throw OomphLibError(error_message_stream.str(),
812 OOMPH_CURRENT_FUNCTION,
813 OOMPH_EXCEPTION_LOCATION);
814 }
815 // Check that the size of x is correct
816 if (x.size() != N)
817 {
818 std::ostringstream error_message_stream;
819 error_message_stream << "The x vector is not the right size. It is "
820 << x.size() << ", it should be " << N << std::endl;
821
822 throw OomphLibError(error_message_stream.str(),
823 OOMPH_CURRENT_FUNCTION,
824 OOMPH_EXCEPTION_LOCATION);
825 }
826#endif
827
828 unsigned long r_n = residual.size();
829 if (r_n != N)
830 {
831 residual.resize(N);
832 }
833
834 // Need to do this in loop over rows
835 for (unsigned i = 0; i < N; i++)
836 {
837 residual[i] = rhs[i];
838 }
839 // Now loop over columns
840 for (unsigned long j = 0; j < N; j++)
841 {
842 for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
843 {
844 unsigned long i = Row_index[k];
845 std::complex<double> a_ij = Value[k];
846 residual[i] -= a_ij * x[j];
847 }
848 }
849 }
850
851 //===================================================================
852 /// Multiply the matrix by the vector x
853 //===================================================================
854 void CCComplexMatrix::multiply(const Vector<std::complex<double>>& x,
855 Vector<std::complex<double>>& soln)
856 {
857#ifdef PARANOID
858 // Check to see if x.size() = ncol()
859 if (x.size() != M)
860 {
861 std::ostringstream error_message_stream;
862 error_message_stream << "The x vector is not the right size. It is "
863 << x.size() << ", it should be " << M << std::endl;
864
865 throw OomphLibError(error_message_stream.str(),
866 OOMPH_CURRENT_FUNCTION,
867 OOMPH_EXCEPTION_LOCATION);
868 }
869#endif
870
871 if (soln.size() != N)
872 {
873 // Resize and initialize the solution vector
874 soln.resize(N);
875 }
876 for (unsigned i = 0; i < N; i++)
877 {
878 soln[i] = 0.0;
879 }
880
881 for (unsigned long j = 0; j < N; j++)
882 {
883 for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
884 {
885 unsigned long i = Row_index[k];
886 std::complex<double> a_ij = Value[k];
887 soln[i] += a_ij * x[j];
888 }
889 }
890 }
891
892
893 //=================================================================
894 /// Multiply the transposed matrix by the vector x: soln=A^T x
895 //=================================================================
897 const Vector<std::complex<double>>& x, Vector<std::complex<double>>& soln)
898 {
899#ifdef PARANOID
900 // Check to see x.size() = nrow()
901 if (x.size() != N)
902 {
903 std::ostringstream error_message_stream;
904 error_message_stream << "The x vector is not the right size. It is "
905 << x.size() << ", it should be " << N << std::endl;
906
907 throw OomphLibError(error_message_stream.str(),
908 OOMPH_CURRENT_FUNCTION,
909 OOMPH_EXCEPTION_LOCATION);
910 }
911#endif
912
913 if (soln.size() != M)
914 {
915 // Resize and initialize the solution vector
916 soln.resize(M);
917 }
918
919 // Initialise the solution
920 for (unsigned long i = 0; i < M; i++)
921 {
922 soln[i] = 0.0;
923 }
924
925 // Matrix vector product
926 for (unsigned long i = 0; i < N; i++)
927 {
928 for (long k = Column_start[i]; k < Column_start[i + 1]; k++)
929 {
930 unsigned long j = Row_index[k];
931 std::complex<double> a_ij = Value[k];
932 soln[j] += a_ij * x[i];
933 }
934 }
935 }
936
937
938 /// /////////////////////////////////////////////////////////////////////
939 /// /////////////////////////////////////////////////////////////////////
940 /// /////////////////////////////////////////////////////////////////////
941
942
943 //===================================================================
944 /// Do LU decomposition and return sign of determinant
945 //===================================================================
947 {
948#ifdef PARANOID
949 if (N != M)
950 {
951 std::ostringstream error_message_stream;
952 error_message_stream << "Can only solve for quadratic matrices\n"
953 << "N, M " << N << " " << M << std::endl;
954
955 throw OomphLibError(error_message_stream.str(),
956 OOMPH_CURRENT_FUNCTION,
957 OOMPH_EXCEPTION_LOCATION);
958 }
959#endif
960
961 // SuperLU expects compressed column format: Set flag for
962 // transpose (0/1) = (false/true)
963 int transpose = 1;
964
965 // Doc (0/1) = (false/true)
966 int doc = 0;
968
969 // Copies so that const-ness is maintained
970 int n_aux = int(N);
971 int nnz_aux = int(Nnz);
972
973 // Integer to hold the sign of the determinant
974 int sign = 0;
975
976 // Just do the lu decompose phase (i=1)
977 int i = 1;
978 sign = superlu_complex(&i,
979 &n_aux,
980 &nnz_aux,
981 0,
982 Value,
984 Row_start,
985 0,
986 &n_aux,
987 &transpose,
988 &doc,
989 &F_factors,
990 &Info);
991 // Return sign of determinant
992 return sign;
993 }
994
995
996 //===================================================================
997 /// Do back-substitution
998 //===================================================================
999 void CRComplexMatrix::lubksub(Vector<std::complex<double>>& rhs)
1000 {
1001#ifdef PARANOID
1002 if (N != rhs.size())
1003 {
1004 std::ostringstream error_message_stream;
1005 error_message_stream << "Size of RHS doesn't match matrix size\n"
1006 << "N, rhs.size() " << N << " " << rhs.size()
1007 << std::endl;
1008
1009 throw OomphLibError(error_message_stream.str(),
1010 OOMPH_CURRENT_FUNCTION,
1011 OOMPH_EXCEPTION_LOCATION);
1012 }
1013 if (N != M)
1014 {
1015 std::ostringstream error_message_stream;
1016 error_message_stream << "Can only solve for quadratic matrices\n"
1017 << "N, M " << N << " " << M << std::endl;
1018
1019 throw OomphLibError(error_message_stream.str(),
1020 OOMPH_CURRENT_FUNCTION,
1021 OOMPH_EXCEPTION_LOCATION);
1022 }
1023#endif
1024
1025 /// RHS vector
1026 std::complex<double>* b = new std::complex<double>[N];
1027
1028 // Copy across
1029 for (unsigned long i = 0; i < N; i++)
1030 {
1031 b[i] = rhs[i];
1032 }
1033
1034 // SuperLU expects compressed column format: Set flag for
1035 // transpose (0/1) = (false/true)
1036 int transpose = 1;
1037
1038 // Doc (0/1) = (false/true)
1039 int doc = 0;
1040 if (Doc_stats_during_solve) doc = 1;
1041
1042 // Number of RHSs
1043 int nrhs = 1;
1044
1045 // Copies so that const-ness is maintained
1046 int n_aux = int(N);
1047 int nnz_aux = int(Nnz);
1048
1049 // SuperLU in three steps (lu decompose, back subst, cleanup)
1050 // Do the solve phase
1051 int i = 2;
1053 &n_aux,
1054 &nnz_aux,
1055 &nrhs,
1056 Value,
1058 Row_start,
1059 b,
1060 &n_aux,
1061 &transpose,
1062 &doc,
1063 &F_factors,
1064 &Info);
1065
1066 // Copy b into rhs vector
1067 for (unsigned long i = 0; i < N; i++)
1068 {
1069 rhs[i] = b[i];
1070 }
1071
1072 // Cleanup
1073 delete[] b;
1074 }
1075
1076
1077 //===================================================================
1078 /// Cleanup memory
1079 //===================================================================
1081 {
1082 // Only bother if we've done an LU solve
1083 if (F_factors != 0)
1084 {
1085 // SuperLU expects compressed column format: Set flag for
1086 // transpose (0/1) = (false/true)
1087 int transpose = 1;
1088
1089 // Doc (0/1) = (false/true)
1090 int doc = 0;
1091 if (Doc_stats_during_solve) doc = 1;
1092
1093 // Copies so that const-ness is maintained
1094 int n_aux = int(N);
1095 int nnz_aux = int(Nnz);
1096
1097 // SuperLU in three steps (lu decompose, back subst, cleanup)
1098 // Flag to indicate which solve step to do (1, 2 or 3)
1099 int i = 3;
1101 &n_aux,
1102 &nnz_aux,
1103 0,
1104 Value,
1106 Row_start,
1107 0,
1108 &n_aux,
1109 &transpose,
1110 &doc,
1111 &F_factors,
1112 &Info);
1113 }
1114 }
1115
1116
1117 //=================================================================
1118 /// Find the residulal to x of Ax=b, ie r=b-Ax
1119 //=================================================================
1120 void CRComplexMatrix::residual(const Vector<std::complex<double>>& x,
1121 const Vector<std::complex<double>>& rhs,
1122 Vector<std::complex<double>>& residual)
1123 {
1124#ifdef PARANOID
1125 // Check that size of rhs = nrow()
1126 if (rhs.size() != N)
1127 {
1128 std::ostringstream error_message_stream;
1129 error_message_stream << "The rhs vector is not the right size. It is "
1130 << rhs.size() << ", it should be " << N << std::endl;
1131
1132 throw OomphLibError(error_message_stream.str(),
1133 OOMPH_CURRENT_FUNCTION,
1134 OOMPH_EXCEPTION_LOCATION);
1135 }
1136 // Check that the size of x is correct
1137 if (x.size() != M)
1138 {
1139 std::ostringstream error_message_stream;
1140 error_message_stream << "The x vector is not the right size. It is "
1141 << x.size() << ", it should be " << M << std::endl;
1142
1143 throw OomphLibError(error_message_stream.str(),
1144 OOMPH_CURRENT_FUNCTION,
1145 OOMPH_EXCEPTION_LOCATION);
1146 }
1147#endif
1148
1149 if (residual.size() != N)
1150 {
1151 residual.resize(N);
1152 }
1153
1154 for (unsigned long i = 0; i < N; i++)
1155 {
1156 residual[i] = rhs[i];
1157 for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
1158 {
1159 unsigned long j = Column_index[k];
1160 std::complex<double> a_ij = Value[k];
1161 residual[i] -= a_ij * x[j];
1162 }
1163 }
1164 }
1165
1166 //=================================================================
1167 /// Multiply the matrix by the vector x
1168 //=================================================================
1169 void CRComplexMatrix::multiply(const Vector<std::complex<double>>& x,
1170 Vector<std::complex<double>>& soln)
1171 {
1172#ifdef PARANOID
1173 // Check to see x.size() = ncol()
1174 if (x.size() != M)
1175 {
1176 std::ostringstream error_message_stream;
1177 error_message_stream << "The x vector is not the right size. It is "
1178 << x.size() << ", it should be " << M << std::endl;
1179
1180 throw OomphLibError(error_message_stream.str(),
1181 OOMPH_CURRENT_FUNCTION,
1182 OOMPH_EXCEPTION_LOCATION);
1183 }
1184#endif
1185
1186 if (soln.size() != N)
1187 {
1188 // Resize and initialize the solution vector
1189 soln.resize(N);
1190 }
1191 for (unsigned long i = 0; i < N; i++)
1192 {
1193 soln[i] = 0.0;
1194 for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
1195 {
1196 unsigned long j = Column_index[k];
1197 std::complex<double> a_ij = Value[k];
1198 soln[i] += a_ij * x[j];
1199 }
1200 }
1201 }
1202
1203
1204 //=================================================================
1205 /// Multiply the transposed matrix by the vector x: soln=A^T x
1206 //=================================================================
1208 const Vector<std::complex<double>>& x, Vector<std::complex<double>>& soln)
1209 {
1210#ifdef PARANOID
1211 // Check to see x.size() = nrow()
1212 if (x.size() != N)
1213 {
1214 std::ostringstream error_message_stream;
1215 error_message_stream << "The x vector is not the right size. It is "
1216 << x.size() << ", it should be " << N << std::endl;
1217
1218 throw OomphLibError(error_message_stream.str(),
1219 OOMPH_CURRENT_FUNCTION,
1220 OOMPH_EXCEPTION_LOCATION);
1221 }
1222#endif
1223
1224 if (soln.size() != M)
1225 {
1226 // Resize and initialize the solution vector
1227 soln.resize(M);
1228 }
1229
1230 // Initialise the solution
1231 for (unsigned long i = 0; i < M; i++)
1232 {
1233 soln[i] = 0.0;
1234 }
1235
1236 // Matrix vector product
1237 for (unsigned long i = 0; i < N; i++)
1238 {
1239 for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
1240 {
1241 unsigned long j = Column_index[k];
1242 std::complex<double> a_ij = Value[k];
1243 soln[j] += a_ij * x[i];
1244 }
1245 }
1246 }
1247} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
int Info
Info flag for the SuperLU solver.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
int ludecompose()
LU decomposition using SuperLU.
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
void clean_up_memory()
LU clean up (perhaps this should happen in the destructor)
bool Doc_stats_during_solve
Flag to indicate if stats are to be displayed during solution of linear system with SuperLU.
void * F_factors
Storage for the LU factors as required by SuperLU.
void lubksub(Vector< std::complex< double > > &rhs)
LU back solve for given RHS.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &b, Vector< std::complex< double > > &residual)
Find the residulal to x of Ax=b, ie r=b-Ax.
int * Column_start
Start index for column.
Definition: matrices.h:2779
void * F_factors
Storage for the LU factors as required by SuperLU.
bool Doc_stats_during_solve
Flag to indicate if stats are to be displayed during solution of linear system with SuperLU.
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
int ludecompose()
LU decomposition using SuperLU.
void lubksub(Vector< std::complex< double > > &rhs)
LU back solve for given RHS.
int Info
Info flag for the SuperLU solver.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &b, Vector< std::complex< double > > &residual)
Find the residual to x of Ax=b, i.e. r=b-Ax.
void clean_up_memory()
LU clean up (perhaps this should happen in the destructor)
int * Row_start
Start index for row.
Definition: matrices.h:873
virtual int ludecompose()
LU decomposition of the matrix using the appropriate linear solver. Return the sign of the determinan...
virtual void solve(Vector< std::complex< double > > &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual void lubksub(Vector< std::complex< double > > &rhs)
LU backsubstitue a previously LU-decomposed matrix; The passed rhs will be over-written with the solu...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
bool Overwrite_matrix_storage
Boolean flag used to decided whether the LU decomposition is stored separately, or not.
virtual ~DenseComplexMatrix()
Destructor, delete the storage for Index vector and LU storage (if any)
Vector< long > * Index
Pointer to storage for the index of permutations in the LU solve.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &rhs, Vector< std::complex< double > > &residual)
Find the residual of Ax=b, ie r=b-Ax for the "solution" x.
int ludecompose()
Overload the LU decomposition. For this storage scheme the matrix storage will be over-writeen by its...
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
void delete_lu_factors()
Function that deletes the storage for the LU_factors, if it is not the same as the matrix storage.
void lubksub(Vector< std::complex< double > > &rhs)
Overload the LU back substitution. Note that the rhs will be overwritten with the solution vector.
std::complex< double > * LU_factors
Pointer to storage for the LU decomposition.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
unsigned long N
Number of rows.
Definition: matrices.h:392
std::complex< double > & entry(const unsigned long &i, const unsigned long &j)
The access function that will be called by the read-write round-bracket operator.
Definition: matrices.h:447
std::complex< double > * Matrixdata
Internal representation of matrix as a pointer to data.
Definition: matrices.h:389
unsigned long M
Number of columns.
Definition: matrices.h:395
An OomphLibError object which should be thrown when an run-time error is encountered....
unsigned long Nnz
Number of non-zero values (i.e. size of Value array)
Definition: matrices.h:574
T * Value
Internal representation of the matrix values, a pointer.
Definition: matrices.h:565
unsigned long N
Number of rows.
Definition: matrices.h:568
unsigned long M
Number of columns.
Definition: matrices.h:571
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
int superlu_complex(int *, int *, int *, int *, std::complex< double > *, int *, int *, std::complex< double > *, int *, int *, int *, void *, int *)
///////////////////////////////////////////////////////////////////// ///////////////////////////////...