00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include "CmplxQR.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034
00035 extern "C"
00036 {
00037 F77_RET_T
00038 F77_FUNC (zgeqrf, ZGEQRF) (const int&, const int&, Complex*,
00039 const int&, Complex*, Complex*,
00040 const int&, int&);
00041
00042 F77_RET_T
00043 F77_FUNC (zungqr, ZUNGQR) (const int&, const int&, const int&,
00044 Complex*, const int&, Complex*,
00045 Complex*, const int&, int&);
00046 }
00047
00048 ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type)
00049 : q (), r ()
00050 {
00051 init (a, qr_type);
00052 }
00053
00054 void
00055 ComplexQR::init (const ComplexMatrix& a, QR::type qr_type)
00056 {
00057 int m = a.rows ();
00058 int n = a.cols ();
00059
00060 if (m == 0 || n == 0)
00061 {
00062 (*current_liboctave_error_handler)
00063 ("ComplexQR must have non-empty matrix");
00064 return;
00065 }
00066
00067 int min_mn = m < n ? m : n;
00068
00069 Array<Complex> tau (min_mn);
00070 Complex *ptau = tau.fortran_vec ();
00071
00072 int lwork = 32*n;
00073 Array<Complex> work (lwork);
00074 Complex *pwork = work.fortran_vec ();
00075
00076 int info = 0;
00077
00078 ComplexMatrix A_fact;
00079 if (m > n && qr_type != QR::economy)
00080 {
00081 A_fact.resize (m, m);
00082 A_fact.insert (a, 0, 0);
00083 }
00084 else
00085 A_fact = a;
00086
00087 Complex *tmp_data = A_fact.fortran_vec ();
00088
00089 F77_XFCN (zgeqrf, ZGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info));
00090
00091 if (f77_exception_encountered)
00092 (*current_liboctave_error_handler) ("unrecoverable error in zgeqrf");
00093 else
00094 {
00095 if (qr_type == QR::raw)
00096 {
00097 for (int j = 0; j < min_mn; j++)
00098 {
00099 int limit = j < min_mn - 1 ? j : min_mn - 1;
00100 for (int i = limit + 1; i < m; i++)
00101 A_fact.elem (i, j) *= tau.elem (j);
00102 }
00103
00104 r = A_fact;
00105
00106 if (m > n)
00107 r.resize (m, n);
00108 }
00109 else
00110 {
00111 int n2 = (qr_type == QR::economy) ? min_mn : m;
00112
00113 if (qr_type == QR::economy && m > n)
00114 r.resize (n, n, 0.0);
00115 else
00116 r.resize (m, n, 0.0);
00117
00118 for (int j = 0; j < n; j++)
00119 {
00120 int limit = j < min_mn-1 ? j : min_mn-1;
00121 for (int i = 0; i <= limit; i++)
00122 r.elem (i, j) = A_fact.elem (i, j);
00123 }
00124
00125 lwork = 32 * n2;
00126 work.resize (lwork);
00127 Complex *pwork2 = work.fortran_vec ();
00128
00129 F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
00130 pwork2, lwork, info));
00131
00132 if (f77_exception_encountered)
00133 (*current_liboctave_error_handler)
00134 ("unrecoverable error in zungqr");
00135 else
00136 {
00137 q = A_fact;
00138 q.resize (m, n2);
00139 }
00140 }
00141 }
00142 }
00143
00144
00145
00146
00147
00148