MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
complex_densemat.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #include "complex_densemat.hpp"
13 #include <complex>
14 
15 #ifdef MFEM_USE_LAPACK
16 extern "C" void
17 zgetrf_(int *, int *, std::complex<double> *, int *, int *, int *);
18 extern "C" void
19 zgetrs_(char *, int *, int *, std::complex<double> *, int *, int *,
20  std::complex<double> *, int *, int *);
21 extern "C" void
22 zgetri_(int *, std::complex<double> *, int *, int *,
23  std::complex<double> *, int *, int *);
24 extern "C" void
25 ztrsm_(char *, char *, char *, char *, int *, int *, std::complex<double> *,
26  std::complex<double> *, int *, std::complex<double> *, int *);
27 extern "C" void
28 zpotrf_(char *, int *, std::complex<double> *, int *, int *);
29 
30 extern "C" void
31 ztrtrs_(char *, char*, char *, int *, int *, std::complex<double> *, int *,
32  std::complex<double> *, int *, int *);
33 extern "C" void
34 zpotri_(char *, int *, std::complex<double> *, int*, int *);
35 
36 extern "C" void
37 zpotrs_(char *, int *, int *, std::complex<double> *, int *,
38  std::complex<double> *, int *, int *);
39 #endif
40 
41 namespace mfem
42 {
43 
45 {
46  MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
47  return dynamic_cast<DenseMatrix &>(*Op_Real_);
48 }
49 
51 {
52  MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
53  return dynamic_cast<DenseMatrix &>(*Op_Imag_);
54 }
55 
57 {
58  MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
59  return dynamic_cast<const DenseMatrix &>(*Op_Real_);
60 }
61 
63 {
64  MFEM_ASSERT(Op_Imag_, "ComplexDenseMatrix has no imaginary part!");
65  return dynamic_cast<const DenseMatrix &>(*Op_Imag_);
66 }
67 
69 {
70  int h = height/2;
71  int w = width/2;
72  DenseMatrix * A = new DenseMatrix(2*h,2*w);
73  double * data = A->Data();
74  double * data_r = nullptr;
75  double * data_i = nullptr;
76 
77  const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
78  *A = 0.;
79  if (hasRealPart())
80  {
81  data_r = real().Data();
82  for (int j = 0; j<w; j++)
83  {
84  for (int i = 0; i<h; i++)
85  {
86  data[i+j*height] = data_r[i+j*h];
87  data[i+h+(j+h)*height] = factor*data_r[i+j*h];
88  }
89  }
90  }
91  if (hasImagPart())
92  {
93  data_i = imag().Data();
94  for (int j = 0; j<w; j++)
95  {
96  for (int i = 0; i<h; i++)
97  {
98  data[i+h+j*height] = factor*data_i[i+j*h];
99  data[i+(j+h)*height] = -data_i[i+j*h];
100  }
101  }
102  }
103  return A;
104 }
105 
107 {
108  MFEM_VERIFY(height == width, "Matrix has to be square");
109 
110  // complex data
111  int h = height/2;
112  int w = width/2;
113  std::complex<double> * data = new std::complex<double>[h*w];
114 
115  // copy data
116  if (hasRealPart() && hasImagPart())
117  {
118  double * data_r = real().Data();
119  double * data_i = imag().Data();
120  for (int i = 0; i < h*w; i++)
121  {
122  data[i] = std::complex<double> (data_r[i], data_i[i]);
123  }
124  }
125  else if (hasRealPart())
126  {
127  double * data_r = real().Data();
128  for (int i = 0; i < h*w; i++)
129  {
130  data[i] = std::complex<double> (data_r[i], 0.);
131  }
132  }
133  else if (hasImagPart())
134  {
135  double * data_i = imag().Data();
136  for (int i = 0; i < h*w; i++)
137  {
138  data[i] = std::complex<double> (0., data_i[i]);
139  }
140  }
141  else
142  {
143  MFEM_ABORT("ComplexDenseMatrix has neither only a real nor an imag part");
144  }
145 
146 #ifdef MFEM_USE_LAPACK
147  int *ipiv = new int[w];
148  int lwork = -1;
149  std::complex<double> qwork, *work;
150  int info;
151 
152  zgetrf_(&w, &w, data, &w, ipiv, &info);
153  if (info)
154  {
155  mfem_error("DenseMatrix::Invert() : Error in ZGETRF");
156  }
157 
158  zgetri_(&w, data, &w, ipiv, &qwork, &lwork, &info);
159 
160  lwork = (int) qwork.real();
161  work = new std::complex<double>[lwork];
162 
163  zgetri_(&w, data, &w, ipiv, work, &lwork, &info);
164 
165  if (info)
166  {
167  mfem_error("DenseMatrix::Invert() : Error in ZGETRI");
168  }
169 
170  delete [] work;
171  delete [] ipiv;
172 
173 #else
174  // compiling without LAPACK
175  int c, i, j, n = w;
176  double a, b;
177  Array<int> piv(n);
178  std::complex<double> ac,bc;
179 
180  for (c = 0; c < n; c++)
181  {
182  a = std::abs(data[c+c*h]);
183  i = c;
184  for (j = c + 1; j < n; j++)
185  {
186  b = std::abs(data[j+c*h]);
187  if (a < b)
188  {
189  a = b;
190  i = j;
191  }
192  }
193  if (a == 0.0)
194  {
195  mfem_error("DenseMatrix::Invert() : singular matrix");
196  }
197  piv[c] = i;
198  for (j = 0; j < n; j++)
199  {
200  mfem::Swap<std::complex<double>>(data[c+j*h], data[i+j*h]);
201  }
202 
203  ac = data[c+c*h] = 1.0 / data[c+c*h];
204  for (j = 0; j < c; j++)
205  {
206  data[c+j*h] *= ac;
207  }
208  for (j++; j < n; j++)
209  {
210  data[c+j*h] *= ac;
211  }
212  for (i = 0; i < c; i++)
213  {
214  data[i+c*h] = ac * (bc = -data[i+c*h]);
215  for (j = 0; j < c; j++)
216  {
217  data[i+j*h] += bc * data[c+j*h];
218  }
219  for (j++; j < n; j++)
220  {
221  data[i+j*h] += bc * data[c+j*h];
222  }
223  }
224  for (i++; i < n; i++)
225  {
226  data[i+c*h] = ac * (bc = -data[i+c*h]);
227  for (j = 0; j < c; j++)
228  {
229  data[i+j*h] += bc * data[c+j*h];
230  }
231  for (j++; j < n; j++)
232  {
233  data[i+j*h] += bc * data[c+j*h];
234  }
235  }
236  }
237 
238  for (c = n - 1; c >= 0; c--)
239  {
240  j = piv[c];
241  for (i = 0; i < n; i++)
242  {
243  mfem::Swap<std::complex<double>>(data[i+c*h], data[i+j*h]);
244  }
245  }
246 
247 #endif
248 
249  DenseMatrix * C_r = new DenseMatrix(h);
250  DenseMatrix * C_i = new DenseMatrix(h);
251 
252  double * datac_r = C_r->Data();
253  double * datac_i = C_i->Data();
254 
255  for (int k = 0; k < h*w; k++)
256  {
257  datac_r[k] = data[k].real();
258  datac_i[k] = data[k].imag();
259  }
260  delete [] data;
261  return new ComplexDenseMatrix(C_r,C_i,true,true);
262 
263 }
264 
266  const ComplexDenseMatrix &B)
267 {
268  // C = C_r + i C_i = (A_r + i * A_i) * (B_r + i * B_i)
269  // = A_r * B_r - A_i B_i + i (A_r * B_i + A_i * B_r)
270 
271  int h = A.Height()/2;
272  int w = B.Width()/2;
273 
274  MFEM_VERIFY(A.Width() == B.Height(), "Incompatible matrix dimensions");
275 
276  //only real case (imag is null)
277  DenseMatrix * C_r = nullptr;
278  DenseMatrix * C_i = nullptr;
279  if ((A.hasRealPart() && B.hasRealPart()) ||
280  (A.hasImagPart() && B.hasImagPart()))
281  {
282  C_r = new DenseMatrix(h,w);
283  }
284  if ((A.hasRealPart() && B.hasImagPart()) ||
285  (A.hasImagPart() && B.hasRealPart()))
286  {
287  C_i = new DenseMatrix(h,w);
288  }
289 
290  MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
291 
292  if (A.hasRealPart() && B.hasRealPart())
293  {
294  Mult(A.real(), B.real(),*C_r);
295  }
296  if (A.hasImagPart() && B.hasImagPart())
297  {
298  if (A.hasRealPart() && B.hasRealPart())
299  {
300  AddMult_a(-1.,A.imag(), B.imag(),*C_r);
301  }
302  else
303  {
304  Mult(A.imag(), B.imag(),*C_r);
305  }
306  }
307 
308  if (A.hasRealPart() && B.hasImagPart())
309  {
310  Mult(A.real(), B.imag(),*C_i);
311  }
312 
313  if (A.hasImagPart() && B.hasRealPart())
314  {
315  if (A.hasRealPart() && B.hasImagPart())
316  {
317  AddMult(A.imag(), B.real(),*C_i);
318  }
319  else
320  {
321  Mult(A.imag(), B.real(),*C_i);
322  }
323  }
324 
325  return new ComplexDenseMatrix(C_r,C_i,true,true);
326 }
327 
329  const ComplexDenseMatrix &B)
330 {
331  // C = C_r + i C_i = (A_r^t - i * A_i^t) * (B_r + i * B_i)
332  // = A_r^t * B_r + A_i^t * B_i + i (A_r^t * B_i - A_i^t * B_r)
333 
334  int h = A.Width()/2;
335  int w = B.Width()/2;
336 
337  MFEM_VERIFY(A.Height() == B.Height(), "Incompatible matrix dimensions");
338 
339  //only real case (imag is null)
340  DenseMatrix * C_r = nullptr;
341  DenseMatrix * C_i = nullptr;
342  if ((A.hasRealPart() && B.hasRealPart()) ||
343  (A.hasImagPart() && B.hasImagPart()))
344  {
345  C_r = new DenseMatrix(h,w);
346  }
347  if ((A.hasRealPart() && B.hasImagPart()) ||
348  (A.hasImagPart() && B.hasRealPart()))
349  {
350  C_i = new DenseMatrix(h,w);
351  }
352 
353  MFEM_VERIFY(C_r || C_i, "Both real and imag parts are null");
354 
355  if (A.hasRealPart() && B.hasRealPart())
356  {
357  MultAtB(A.real(), B.real(),*C_r);
358  }
359  if (A.hasImagPart() && B.hasImagPart())
360  {
361  if (A.hasRealPart() && B.hasRealPart())
362  {
363  DenseMatrix tempC_r(h,w);
364  MultAtB(A.imag(), B.imag(),tempC_r);
365  (*C_r) += tempC_r;
366  }
367  else
368  {
369  MultAtB(A.imag(), B.imag(),*C_r);
370  }
371  }
372 
373  if (A.hasRealPart() && B.hasImagPart())
374  {
375  MultAtB(A.real(), B.imag(),*C_i);
376  }
377 
378  if (A.hasImagPart() && B.hasRealPart())
379  {
380  if (A.hasRealPart() && B.hasImagPart())
381  {
382  DenseMatrix tempC_i(h,w);
383  MultAtB(A.imag(), B.real(),tempC_i);
384  (*C_i) -= tempC_i;
385  }
386  else
387  {
388  MultAtB(A.imag(), B.real(),*C_i);
389  }
390  }
391 
392  return new ComplexDenseMatrix(C_r,C_i,true,true);
393 }
394 
395 std::complex<double> * ComplexFactors::RealToComplex
396 (int m, const double * x_r, const double * x_i) const
397 {
398  std::complex<double> * x = new std::complex<double>[m];
399  if (x_r && x_i)
400  {
401  for (int i = 0; i<m; i++)
402  {
403  x[i] = std::complex<double>(x_r[i], x_i[i]);
404  }
405  }
406  else if (data_r)
407  {
408  for (int i = 0; i<m; i++)
409  {
410  x[i] = std::complex<double>(x_r[i], 0.);
411  }
412  }
413  else if (data_i)
414  {
415  for (int i = 0; i<m; i++)
416  {
417  x[i] = std::complex<double>(0., x_i[i]);
418  }
419  }
420  else
421  {
422  MFEM_ABORT("ComplexFactors::RealToComplex:both real and imag part are null");
423  return nullptr;
424  }
425  return x;
426 }
427 
428 void ComplexFactors::ComplexToReal(int m, const std::complex<double> * x,
429  double * x_r, double * x_i) const
430 {
431  for (int i = 0; i<m; i++)
432  {
433  x_r[i] = x[i].real();
434  x_i[i] = x[i].imag();
435  }
436 }
437 
438 
440 {
441  if (data) { return; }
442  MFEM_VERIFY(data_r || data_i, "ComplexFactors data not set");
444 }
445 
446 
447 bool ComplexLUFactors::Factor(int m, double TOL)
448 {
449  SetComplexData(m*m);
450 #ifdef MFEM_USE_LAPACK
451  int info = 0;
452  MFEM_VERIFY(data, "Matrix data not set");
453  if (m) { zgetrf_(&m, &m, data, &m, ipiv, &info); }
454  return info == 0;
455 #else
456  // compiling without LAPACK
457  std::complex<double> *data_ptr = this->data;
458  for (int i = 0; i < m; i++)
459  {
460  // pivoting
461  {
462  int piv = i;
463  double a = std::abs(data_ptr[piv+i*m]);
464  for (int j = i+1; j < m; j++)
465  {
466  const double b = std::abs(data_ptr[j+i*m]);
467  if (b > a)
468  {
469  a = b;
470  piv = j;
471  }
472  }
473  ipiv[i] = piv;
474  if (piv != i)
475  {
476  // swap rows i and piv in both L and U parts
477  for (int j = 0; j < m; j++)
478  {
479  mfem::Swap<std::complex<double>>(data_ptr[i+j*m], data_ptr[piv+j*m]);
480  }
481  }
482  }
483 
484  if (abs(data_ptr[i + i*m]) <= TOL)
485  {
486  return false; // failed
487  }
488 
489  const std::complex<double> a_ii_inv = 1.0 / data_ptr[i+i*m];
490  for (int j = i+1; j < m; j++)
491  {
492  data_ptr[j+i*m] *= a_ii_inv;
493  }
494  for (int k = i+1; k < m; k++)
495  {
496  const std::complex<double> a_ik = data_ptr[i+k*m];
497  for (int j = i+1; j < m; j++)
498  {
499  data_ptr[j+k*m] -= a_ik * data_ptr[j+i*m];
500  }
501  }
502  }
503 #endif
504 
505  return true; // success
506 }
507 
508 std::complex<double> ComplexLUFactors::Det(int m) const
509 {
510  std::complex<double> det(1.0,0.);
511  for (int i=0; i<m; i++)
512  {
513  if (ipiv[i] != i-ipiv_base)
514  {
515  det *= -data[m * i + i];
516  }
517  else
518  {
519  det *= data[m * i + i];
520  }
521  }
522  return det;
523 }
524 
525 void ComplexLUFactors::Mult(int m, int n, double *X_r, double * X_i) const
526 {
527  std::complex<double> * x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
528  for (int k = 0; k < n; k++)
529  {
530  // X <- U X
531  for (int i = 0; i < m; i++)
532  {
533  std::complex<double> x_i = x[i] * data[i+i*m];
534  for (int j = i+1; j < m; j++)
535  {
536  x_i += x[j] * data[i+j*m];
537  }
538  x[i] = x_i;
539  }
540  // X <- L X
541  for (int i = m-1; i >= 0; i--)
542  {
543  std::complex<double> x_i = x[i];
544  for (int j = 0; j < i; j++)
545  {
546  x_i += x[j] * data[i+j*m];
547  }
548  x[i] = x_i;
549  }
550  // X <- P^{-1} X
551  for (int i = m-1; i >= 0; i--)
552  {
553  mfem::Swap<std::complex<double>>(x[i], x[ipiv[i]-ipiv_base]);
554  }
555  x += m;
556  }
557  ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
558  delete [] x;
559 }
560 
561 void ComplexLUFactors::LSolve(int m, int n, double *X_r, double * X_i) const
562 {
563  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
564  std::complex<double> *x = X;
565  for (int k = 0; k < n; k++)
566  {
567  // X <- P X
568  for (int i = 0; i < m; i++)
569  {
570  mfem::Swap<std::complex<double>>(x[i], x[ipiv[i]-ipiv_base]);
571  }
572  // X <- L^{-1} X
573  for (int j = 0; j < m; j++)
574  {
575  const std::complex<double> x_j = x[j];
576  for (int i = j+1; i < m; i++)
577  {
578  x[i] -= data[i+j*m] * x_j;
579  }
580  }
581  x += m;
582  }
583  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
584  delete [] X;
585 }
586 
587 void ComplexLUFactors::USolve(int m, int n, double *X_r, double * X_i) const
588 {
589  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
590  std::complex<double> *x = X;
591  // X <- U^{-1} X
592  for (int k = 0; k < n; k++)
593  {
594  for (int j = m-1; j >= 0; j--)
595  {
596  const std::complex<double> x_j = ( x[j] /= data[j+j*m] );
597  for (int i = 0; i < j; i++)
598  {
599  x[i] -= data[i+j*m] * x_j;
600  }
601  }
602  x += m;
603  }
604  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
605  delete [] X;
606 }
607 
608 void ComplexLUFactors::Solve(int m, int n, double *X_r, double * X_i) const
609 {
610 #ifdef MFEM_USE_LAPACK
611  std::complex<double> * x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
612  char trans = 'N';
613  int info = 0;
614  if (m > 0 && n > 0) { zgetrs_(&trans, &m, &n, data, &m, ipiv, x, &m, &info); }
615  MFEM_VERIFY(!info, "LAPACK: error in ZGETRS");
616  ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
617  delete [] x;
618 #else
619  // compiling without LAPACK
620  LSolve(m, n, X_r, X_i);
621  USolve(m, n, X_r, X_i);
622 #endif
623 }
624 
625 void ComplexLUFactors::RightSolve(int m, int n, double *X_r, double * X_i) const
626 {
627  std::complex<double> * X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
628  std::complex<double> * x = X;
629 #ifdef MFEM_USE_LAPACK
630  char n_ch = 'N', side = 'R', u_ch = 'U', l_ch = 'L';
631  if (m > 0 && n > 0)
632  {
633  std::complex<double> alpha(1.0,0.0);
634  ztrsm_(&side,&u_ch,&n_ch,&n_ch,&n,&m,&alpha,data,&m,X,&n);
635  ztrsm_(&side,&l_ch,&n_ch,&u_ch,&n,&m,&alpha,data,&m,X,&n);
636  }
637 #else
638  // compiling without LAPACK
639  // X <- X U^{-1}
640  for (int k = 0; k < n; k++)
641  {
642  for (int j = 0; j < m; j++)
643  {
644  const std::complex<double> x_j = ( x[j*n] /= data[j+j*m]);
645  for (int i = j+1; i < m; i++)
646  {
647  x[i*n] -= data[j + i*m] * x_j;
648  }
649  }
650  ++x;
651  }
652 
653  // X <- X L^{-1}
654  x = X;
655  for (int k = 0; k < n; k++)
656  {
657  for (int j = m-1; j >= 0; j--)
658  {
659  const std::complex<double> x_j = x[j*n];
660  for (int i = 0; i < j; i++)
661  {
662  x[i*n] -= data[j + i*m] * x_j;
663  }
664  }
665  ++x;
666  }
667 #endif
668  // X <- X P
669  x = X;
670  for (int k = 0; k < n; k++)
671  {
672  for (int i = m-1; i >= 0; --i)
673  {
674  mfem::Swap<std::complex<double>>(x[i*n], x[(ipiv[i]-ipiv_base)*n]);
675  }
676  ++x;
677  }
678  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
679  delete [] X;
680 }
681 
682 void ComplexLUFactors::GetInverseMatrix(int m, double *X_r, double *X_i) const
683 {
684  // A^{-1} = U^{-1} L^{-1} P
685  // X <- U^{-1} (set only the upper triangular part of X)
686  std::complex<double> * X = ComplexFactors::RealToComplex(m*m,X_r,X_i);
687  std::complex<double> * x = X;
688  for (int k = 0; k < m; k++)
689  {
690  const std::complex<double> minus_x_k = -( x[k] = 1.0/data[k+k*m] );
691  for (int i = 0; i < k; i++)
692  {
693  x[i] = data[i+k*m] * minus_x_k;
694  }
695  for (int j = k-1; j >= 0; j--)
696  {
697  const std::complex<double> x_j = ( x[j] /= data[j+j*m] );
698  for (int i = 0; i < j; i++)
699  {
700  x[i] -= data[i+j*m] * x_j;
701  }
702  }
703  x += m;
704  }
705  // X <- X L^{-1} (use input only from the upper triangular part of X)
706  {
707  int k = m-1;
708  for (int j = 0; j < k; j++)
709  {
710  const std::complex<double> minus_L_kj = -data[k+j*m];
711  for (int i = 0; i <= j; i++)
712  {
713  X[i+j*m] += X[i+k*m] * minus_L_kj;
714  }
715  for (int i = j+1; i < m; i++)
716  {
717  X[i+j*m] = X[i+k*m] * minus_L_kj;
718  }
719  }
720  }
721  for (int k = m-2; k >= 0; k--)
722  {
723  for (int j = 0; j < k; j++)
724  {
725  const std::complex<double> L_kj = data[k+j*m];
726  for (int i = 0; i < m; i++)
727  {
728  X[i+j*m] -= X[i+k*m] * L_kj;
729  }
730  }
731  }
732  // X <- X P
733  for (int k = m-1; k >= 0; k--)
734  {
735  const int piv_k = ipiv[k]-ipiv_base;
736  if (k != piv_k)
737  {
738  for (int i = 0; i < m; i++)
739  {
740  Swap<std::complex<double>>(X[i+k*m], X[i+piv_k*m]);
741  }
742  }
743  }
744  ComplexFactors::ComplexToReal(m*m,X,X_r,X_i);
745  delete [] X;
746 }
747 
748 
749 bool ComplexCholeskyFactors::Factor(int m, double TOL)
750 {
751  SetComplexData(m*m);
752 #ifdef MFEM_USE_LAPACK
753  int info = 0;
754  char uplo = 'L';
755  MFEM_VERIFY(data, "Matrix data not set");
756  if (m) {zpotrf_(&uplo, &m, data, &m, &info);}
757  return info == 0;
758 #else
759  // Cholesky–Crout algorithm
760  for (int j = 0; j<m; j++)
761  {
762  std::complex<double> a(0.,0.);
763  for (int k = 0; k<j; k++)
764  {
765  a+=data[j+k*m]*std::conj(data[j+k*m]);
766  }
767 
768  MFEM_VERIFY((data[j+j*m] - a).real() > 0.,
769  "CholeskyFactors::Factor: The matrix is not SPD");
770 
771  data[j+j*m] = std::sqrt((data[j+j*m] - a).real());
772 
773  if (data[j + j*m].real() <= TOL) { return false; }
774 
775  for (int i = j+1; i<m; i++)
776  {
777  a = std::complex<double>(0.,0.);
778  for (int k = 0; k<j; k++)
779  {
780  a+= data[i+k*m]*std::conj(data[j+k*m]);
781  }
782  data[i+j*m] = 1./data[j+m*j]*(data[i+j*m] - a);
783  }
784  }
785  return true; // success
786 #endif
787 }
788 
789 std::complex<double> ComplexCholeskyFactors::Det(int m) const
790 {
791  std::complex<double> det(1.0,0.0);
792  for (int i=0; i<m; i++)
793  {
794  det *= data[i + i*m];
795  }
796  return det;
797 }
798 
799 void ComplexCholeskyFactors::LMult(int m, int n, double * X_r,
800  double * X_i) const
801 {
802  // X <- L X
803  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
804  std::complex<double> *x = X;
805  for (int k = 0; k < n; k++)
806  {
807  for (int j = m-1; j >= 0; j--)
808  {
809  std::complex<double> x_j = x[j] * data[j+j*m];
810  for (int i = 0; i < j; i++)
811  {
812  x_j += x[i] * data[j+i*m];
813  }
814  x[j] = x_j;
815  }
816  x += m;
817  }
818  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
819  delete [] X;
820 }
821 
822 void ComplexCholeskyFactors::UMult(int m, int n, double * X_r,
823  double * X_i) const
824 {
825  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
826  std::complex<double> *x = X;
827  for (int k = 0; k < n; k++)
828  {
829  for (int i = 0; i < m; i++)
830  {
831  std::complex<double> x_i = x[i] * data[i+i*m];
832  for (int j = i+1; j < m; j++)
833  {
834  x_i += x[j] * std::conj(data[j+i*m]);
835  }
836  x[i] = x_i;
837  }
838  x += m;
839  }
840  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
841  delete [] X;
842 }
843 
844 void ComplexCholeskyFactors::LSolve(int m, int n, double * X_r,
845  double * X_i) const
846 {
847  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
848  std::complex<double> *x = X;
849 #ifdef MFEM_USE_LAPACK
850  char uplo = 'L';
851  char trans = 'N';
852  char diag = 'N';
853  int info = 0;
854 
855  ztrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
856  MFEM_VERIFY(!info, "ComplexCholeskyFactors:LSolve:: info");
857 #else
858  for (int k = 0; k < n; k++)
859  {
860  // X <- L^{-1} X
861  for (int j = 0; j < m; j++)
862  {
863  const std::complex<double> x_j = (x[j] /= data[j+j*m]);
864  for (int i = j+1; i < m; i++)
865  {
866  x[i] -= data[i+j*m] * x_j;
867  }
868  }
869  x += m;
870  }
871 #endif
872  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
873  delete [] X;
874 }
875 
876 void ComplexCholeskyFactors::USolve(int m, int n, double * X_r,
877  double * X_i) const
878 {
879  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
880  std::complex<double> *x = X;
881 #ifdef MFEM_USE_LAPACK
882 
883  char uplo = 'L';
884  char trans = 'C';
885  char diag = 'N';
886  int info = 0;
887 
888  ztrtrs_(&uplo, &trans, &diag, &m, &n, data, &m, x, &m, &info);
889  MFEM_VERIFY(!info, "ComplexCholeskyFactors:USolve:: info");
890 #else
891  // X <- L^{-t} X
892  for (int k = 0; k < n; k++)
893  {
894  for (int j = m-1; j >= 0; j--)
895  {
896  const std::complex<double> x_j = ( x[j] /= data[j+j*m] );
897  for (int i = 0; i < j; i++)
898  {
899  x[i] -= std::conj(data[j+i*m]) * x_j;
900  }
901  }
902  x += m;
903  }
904 #endif
905  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
906  delete [] X;
907 }
908 
909 void ComplexCholeskyFactors::Solve(int m, int n, double * X_r,
910  double * X_i) const
911 {
912 #ifdef MFEM_USE_LAPACK
913  char uplo = 'L';
914  int info = 0;
915  std::complex<double> *x = ComplexFactors::RealToComplex(m*n,X_r,X_i);
916  zpotrs_(&uplo, &m, &n, data, &m, x, &m, &info);
917  MFEM_VERIFY(!info, "ComplexCholeskyFactors:Solve:: info");
918  ComplexFactors::ComplexToReal(m*n,x,X_r,X_i);
919  delete x;
920 #else
921  LSolve(m, n, X_r,X_i);
922  USolve(m, n, X_r,X_i);
923 #endif
924 }
925 
926 void ComplexCholeskyFactors::RightSolve(int m, int n, double * X_r,
927  double * X_i) const
928 {
929 
930  std::complex<double> *X = ComplexFactors::RealToComplex(m*n,X_r,X_i);
931  std::complex<double> *x = X;
932 #ifdef MFEM_USE_LAPACK
933  char side = 'R';
934  char uplo = 'L';
935  char transt = 'C';
936  char trans = 'N';
937  char diag = 'N';
938 
939  std::complex<double> alpha(1.0,0.0);
940  if (m > 0 && n > 0)
941  {
942  ztrsm_(&side,&uplo,&transt,&diag,&n,&m,&alpha,data,&m,x,&n);
943  ztrsm_(&side,&uplo,&trans,&diag,&n,&m,&alpha,data,&m,x,&n);
944  }
945 #else
946  // X <- X L^{-H}
947  for (int k = 0; k < n; k++)
948  {
949  for (int j = 0; j < m; j++)
950  {
951  const std::complex<double> x_j = ( x[j*n] /= data[j+j*m]);
952  for (int i = j+1; i < m; i++)
953  {
954  x[i*n] -= std::conj(data[i + j*m]) * x_j;
955  }
956  }
957  ++x;
958  }
959  // X <- X L^{-1}
960  x = X;
961  for (int k = 0; k < n; k++)
962  {
963  for (int j = m-1; j >= 0; j--)
964  {
965  const std::complex<double> x_j = (x[j*n] /= data[j+j*m]);
966  for (int i = 0; i < j; i++)
967  {
968  x[i*n] -= data[j + i*m] * x_j;
969  }
970  }
971  ++x;
972  }
973 #endif
974  ComplexFactors::ComplexToReal(m*n,X,X_r,X_i);
975  delete [] X;
976 }
977 
979  double * X_i) const
980 {
981  // A^{-1} = L^{-t} L^{-1}
982  std::complex<double> * X = new std::complex<double>[m*m];
983 #ifdef MFEM_USE_LAPACK
984  // copy the lower triangular part of L to X
985  for (int i = 0; i<m; i++)
986  {
987  for (int j = i; j<m; j++)
988  {
989  X[j+i*m] = data[j+i*m];
990  }
991  }
992  char uplo = 'L';
993  int info = 0;
994  zpotri_(&uplo, &m, X, &m, &info);
995  MFEM_VERIFY(!info, "ComplexCholeskyFactors:GetInverseMatrix:: info");
996  // fill in the upper triangular part
997  for (int i = 0; i<m; i++)
998  {
999  for (int j = i+1; j<m; j++)
1000  {
1001  X[i+j*m] = std::conj(X[j+i*m]);
1002  }
1003  }
1004 #else
1005  // L^-t * L^-1 (in place)
1006  for (int k = 0; k<m; k++)
1007  {
1008  X[k+k*m] = 1./data[k+k*m];
1009  for (int i = k+1; i < m; i++)
1010  {
1011  std::complex<double> s(0.,0.);
1012  for (int j=k; j<i; j++)
1013  {
1014  s -= data[i+j*m] * X[j+k*m]/data[i+i*m];
1015  }
1016  X[i+k*m] = s;
1017  }
1018  }
1019  for (int i = 0; i < m; i++)
1020  {
1021  for (int j = i; j < m; j++)
1022  {
1023  std::complex<double> s(0.,0.);
1024  for (int k=j; k<m; k++)
1025  {
1026  s += X[k+i*m] * std::conj(X[k+j*m]);
1027  }
1028  X[j+i*m] = s;
1029  X[i+j*m] = std::conj(s);
1030  }
1031  }
1032 #endif
1033  ComplexFactors::ComplexToReal(m*m,X,X_r,X_i);
1034  delete [] X;
1035 }
1036 
1037 } // mfem namespace
ComplexDenseMatrix * ComputeInverse()
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
void LMult(int m, int n, double *X_r, double *X_i) const
virtual DenseMatrix & real()
Real or imaginary part accessor methods.
void LSolve(int m, int n, double *X_r, double *X_i) const
virtual bool Factor(int m, double TOL=0.0)
Compute the Cholesky factorization of the current matrix.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
std::complex< double > * RealToComplex(int m, const double *x_r, const double *x_i) const
void UMult(int m, int n, double *X_r, double *X_i) const
std::complex< double > * data
void RightSolve(int m, int n, double *X_r, double *X_i) const
void RightSolve(int m, int n, double *X_r, double *X_i) const
bool hasRealPart() const
Check for existence of real or imaginary part of the operator.
void USolve(int m, int n, double *X_r, double *X_i) const
void ztrtrs_(char *, char *, char *, int *, int *, std::complex< double > *, int *, std::complex< double > *, int *, int *)
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2408
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
void Mult(int m, int n, double *X_r, double *X_i) const
void ztrsm_(char *, char *, char *, char *, int *, int *, std::complex< double > *, std::complex< double > *, int *, std::complex< double > *, int *)
DenseMatrix * GetSystemMatrix() const
void LSolve(int m, int n, double *X_r, double *X_i) const
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void zgetrs_(char *, int *, int *, std::complex< double > *, int *, int *, std::complex< double > *, int *, int *)
ComplexDenseMatrix(DenseMatrix *A_Real, DenseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
virtual void GetInverseMatrix(int m, double *X_r, double *X_i) const
Assuming L.L^H = A factored data of size (m x m), compute X &lt;- A^{-1}.
void zgetrf_(int *, int *, std::complex< double > *, int *, int *, int *)
virtual void Solve(int m, int n, double *X_r, double *X_i) const
virtual void GetInverseMatrix(int m, double *X_r, double *X_i) const
Assuming L.U = P.A factored data of size (m x m), compute X &lt;- A^{-1}.
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual std::complex< double > Det(int m) const
void ComplexToReal(int m, const std::complex< double > *x, double *x_r, double *x_i) const
double a
Definition: lissajous.cpp:41
void USolve(int m, int n, double *X_r, double *X_i) const
static const int ipiv_base
virtual void Solve(int m, int n, double *X_r, double *X_i) const
Native convention for Hermitian operators.
virtual std::complex< double > Det(int m) const
const double alpha
Definition: ex15.cpp:369
void zpotri_(char *, int *, std::complex< double > *, int *, int *)
void zpotrs_(char *, int *, int *, std::complex< double > *, int *, std::complex< double > *, int *, int *)
RefCoord s[3]
void zgetri_(int *, std::complex< double > *, int *, int *, std::complex< double > *, int *, int *)
void zpotrf_(char *, int *, std::complex< double > *, int *, int *)
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2375
virtual DenseMatrix & imag()
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
Specialization of the ComplexOperator built from a pair of Dense Matrices. The purpose of this specia...