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