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