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