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 <cassert>
00032
00033 #include "CmplxQRP.h"
00034 #include "f77-fcn.h"
00035 #include "lo-error.h"
00036
00037 extern "C"
00038 {
00039 F77_RET_T
00040 F77_FUNC (zgeqpf, ZGEQPF) (const int&, const int&, Complex*,
00041 const int&, int*, Complex*, Complex*,
00042 double*, int&);
00043
00044 F77_RET_T
00045 F77_FUNC (zungqr, ZUNGQR) (const int&, const int&, const int&,
00046 Complex*, const int&, Complex*,
00047 Complex*, const int&, int&);
00048 }
00049
00050
00051
00052 ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type)
00053 : ComplexQR (), p ()
00054 {
00055 init (a, qr_type);
00056 }
00057
00058 void
00059 ComplexQRP::init (const ComplexMatrix& a, QR::type qr_type)
00060 {
00061 assert (qr_type != QR::raw);
00062
00063 int m = a.rows ();
00064 int n = a.cols ();
00065
00066 if (m == 0 || n == 0)
00067 {
00068 (*current_liboctave_error_handler)
00069 ("ComplexQR must have non-empty matrix");
00070 return;
00071 }
00072
00073 int min_mn = m < n ? m : n;
00074 Array<Complex> tau (min_mn);
00075 Complex *ptau = tau.fortran_vec ();
00076
00077 int lwork = 3*n > 32*m ? 3*n : 32*m;
00078 Array<Complex> work (lwork);
00079 Complex *pwork = work.fortran_vec ();
00080
00081 int info = 0;
00082
00083 ComplexMatrix A_fact = a;
00084 if (m > n && qr_type != QR::economy)
00085 A_fact.resize (m, m, 0.0);
00086
00087 Complex *tmp_data = A_fact.fortran_vec ();
00088
00089 Array<double> rwork (2*n);
00090 double *prwork = rwork.fortran_vec ();
00091
00092 Array<int> jpvt (n, 0);
00093 int *pjpvt = jpvt.fortran_vec ();
00094
00095
00096
00097 F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork,
00098 prwork, info));
00099
00100 if (f77_exception_encountered)
00101 (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf");
00102 else
00103 {
00104
00105
00106
00107 if (qr_type == QR::economy)
00108 {
00109 p.resize (1, n, 0.0);
00110 for (int j = 0; j < n; j++)
00111 p.elem (0, j) = jpvt.elem (j);
00112 }
00113 else
00114 {
00115 p.resize (n, n, 0.0);
00116 for (int j = 0; j < n; j++)
00117 p.elem (jpvt.elem (j) - 1, j) = 1.0;
00118 }
00119
00120 int n2 = (qr_type == QR::economy) ? min_mn : m;
00121
00122 if (qr_type == QR::economy && m > n)
00123 r.resize (n, n, 0.0);
00124 else
00125 r.resize (m, n, 0.0);
00126
00127 for (int j = 0; j < n; j++)
00128 {
00129 int limit = j < min_mn-1 ? j : min_mn-1;
00130 for (int i = 0; i <= limit; i++)
00131 r.elem (i, j) = A_fact.elem (i, j);
00132 }
00133
00134 F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
00135 pwork, lwork, info));
00136
00137 if (f77_exception_encountered)
00138 (*current_liboctave_error_handler) ("unrecoverable error in zungqr");
00139 else
00140 {
00141 q = A_fact;
00142 q.resize (m, n2);
00143 }
00144 }
00145 }
00146
00147
00148
00149
00150
00151