Mercurial > octave

comparison liboctave/dbleQR.cc @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents e35b034d3523
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
29 #include "lo-error.h" 29 #include "lo-error.h"
30 30
31 extern "C" 31 extern "C"
32 { 32 {
33 F77_RET_T 33 F77_RET_T
34 F77_FUNC (dgeqrf, DGEQRF) (const int&, const int&, double*, const int&, 34 F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&,
35 double*, double*, const int&, int&); 35 double*, double*, const octave_idx_type&, octave_idx_type&);
36 36
37 F77_RET_T 37 F77_RET_T
38 F77_FUNC (dorgqr, DORGQR) (const int&, const int&, const int&, double*, 38 F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*,
39 const int&, double*, double*, const int&, int&); 39 const octave_idx_type&, double*, double*, const octave_idx_type&, octave_idx_type&);
40 } 40 }
41 41
42 QR::QR (const Matrix& a, QR::type qr_type) 42 QR::QR (const Matrix& a, QR::type qr_type)
43 : q (), r () 43 : q (), r ()
44 { 44 {
46 } 46 }
47 47
48 void 48 void
49 QR::init (const Matrix& a, QR::type qr_type) 49 QR::init (const Matrix& a, QR::type qr_type)
50 { 50 {
51 int m = a.rows (); 51 octave_idx_type m = a.rows ();
52 int n = a.cols (); 52 octave_idx_type n = a.cols ();
53 53
54 if (m == 0 || n == 0) 54 if (m == 0 || n == 0)
55 { 55 {
56 (*current_liboctave_error_handler) ("QR must have non-empty matrix"); 56 (*current_liboctave_error_handler) ("QR must have non-empty matrix");
57 return; 57 return;
58 } 58 }
59 59
60 int min_mn = m < n ? m : n; 60 octave_idx_type min_mn = m < n ? m : n;
61 Array<double> tau (min_mn); 61 Array<double> tau (min_mn);
62 double *ptau = tau.fortran_vec (); 62 double *ptau = tau.fortran_vec ();
63 63
64 int lwork = 32*n; 64 octave_idx_type lwork = 32*n;
65 Array<double> work (lwork); 65 Array<double> work (lwork);
66 double *pwork = work.fortran_vec (); 66 double *pwork = work.fortran_vec ();
67 67
68 int info = 0; 68 octave_idx_type info = 0;
69 69
70 Matrix A_fact = a; 70 Matrix A_fact = a;
71 if (m > n && qr_type != QR::economy) 71 if (m > n && qr_type != QR::economy)
72 A_fact.resize (m, m, 0.0); 72 A_fact.resize (m, m, 0.0);
73 73
79 (*current_liboctave_error_handler) ("unrecoverable error in dgeqrf"); 79 (*current_liboctave_error_handler) ("unrecoverable error in dgeqrf");
80 else 80 else
81 { 81 {
82 if (qr_type == QR::raw) 82 if (qr_type == QR::raw)
83 { 83 {
84 for (int j = 0; j < min_mn; j++) 84 for (octave_idx_type j = 0; j < min_mn; j++)
85 { 85 {
86 int limit = j < min_mn - 1 ? j : min_mn - 1; 86 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
87 for (int i = limit + 1; i < m; i++) 87 for (octave_idx_type i = limit + 1; i < m; i++)
88 A_fact.elem (i, j) *= tau.elem (j); 88 A_fact.elem (i, j) *= tau.elem (j);
89 } 89 }
90 90
91 r = A_fact; 91 r = A_fact;
92 92
93 if (m > n) 93 if (m > n)
94 r.resize (m, n); 94 r.resize (m, n);
95 } 95 }
96 else 96 else
97 { 97 {
98 int n2 = (qr_type == QR::economy) ? min_mn : m; 98 octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
99 99
100 if (qr_type == QR::economy && m > n) 100 if (qr_type == QR::economy && m > n)
101 r.resize (n, n, 0.0); 101 r.resize (n, n, 0.0);
102 else 102 else
103 r.resize (m, n, 0.0); 103 r.resize (m, n, 0.0);
104 104
105 for (int j = 0; j < n; j++) 105 for (octave_idx_type j = 0; j < n; j++)
106 { 106 {
107 int limit = j < min_mn-1 ? j : min_mn-1; 107 octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
108 for (int i = 0; i <= limit; i++) 108 for (octave_idx_type i = 0; i <= limit; i++)
109 r.elem (i, j) = tmp_data[m*j+i]; 109 r.elem (i, j) = tmp_data[m*j+i];
110 } 110 }
111 111
112 lwork = 32 * n2; 112 lwork = 32 * n2;
113 work.resize (lwork); 113 work.resize (lwork);