Mercurial > octave

comparison liboctave/dbleQR.cc @ 539:5ec10a984241

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-07-21 22:40:04 by jwe]
author jwe
date Thu, 21 Jul 1994 22:42:09 +0000
parents 3d4b4f0fa5ba
children 714fd17fca28
comparison
equal deleted inserted replaced
538:8e134d3b21c9 539:5ec10a984241
39 int F77_FCN (dgeqrf) (const int*, const int*, double*, const int*, 39 int F77_FCN (dgeqrf) (const int*, const int*, double*, const int*,
40 double*, double*, const int*, int*); 40 double*, double*, const int*, int*);
41 41
42 int F77_FCN (dorgqr) (const int*, const int*, const int*, double*, 42 int F77_FCN (dorgqr) (const int*, const int*, const int*, double*,
43 const int*, double*, double*, const int*, int*); 43 const int*, double*, double*, const int*, int*);
44
45 int F77_FCN (dgeqpf) (const int*, const int*, double*, const int*,
46 int*, double*, double*, int*);
44 } 47 }
45 48
46 QR::QR (const Matrix& a) 49 QR::QR (const Matrix& a, QR::type qr_type)
47 { 50 {
48 int m = a.rows (); 51 int m = a.rows ();
49 int n = a.cols (); 52 int n = a.cols ();
50 53
51 if (m == 0 || n == 0) 54 if (m == 0 || n == 0)
71 74
72 F77_FCN (dgeqrf) (&m, &n, tmp_data, &m, tau, work, &lwork, &info); 75 F77_FCN (dgeqrf) (&m, &n, tmp_data, &m, tau, work, &lwork, &info);
73 76
74 delete [] work; 77 delete [] work;
75 78
76 r.resize (m, n, 0.0); 79 if (qr_type == QR::raw)
77 for (int j = 0; j < n; j++)
78 { 80 {
79 int limit = j < min_mn-1 ? j : min_mn-1; 81 for (int j = 0; j < min_mn; j++)
80 for (int i = 0; i <= limit; i++) 82 {
81 r.elem (i, j) = tmp_data[m*j+i]; 83 int limit = j < min_mn - 1 ? j : min_mn - 1;
84 for (int i = limit + 1; i < m; i++)
85 tmp_data[m*j+i] *= tau[j];
86 }
82 } 87 }
88 else
89 {
90 int m2;
91 if (qr_type == QR::economy && m > n)
92 {
93 m2 = n;
94 r.resize (n, n, 0.0);
95 }
96 else
97 {
98 m2 = m;
99 r.resize (m, n, 0.0);
100 }
83 101
84 lwork = 32*m; 102 for (int j = 0; j < n; j++)
85 work = new double[lwork]; 103 {
104 int limit = j < min_mn-1 ? j : min_mn-1;
105 for (int i = 0; i <= limit; i++)
106 r.elem (i, j) = tmp_data[m*j+i];
107 }
86 108
87 F77_FCN (dorgqr) (&m, &m, &min_mn, tmp_data, &m, tau, work, &lwork, &info); 109 lwork = 32*m;
110 work = new double[lwork];
88 111
89 q = Matrix (tmp_data, m, m); 112 F77_FCN (dorgqr) (&m, &m, &min_mn, tmp_data, &m, tau, work,
113 &lwork, &info);
90 114
91 delete [] tau; 115 q = Matrix (tmp_data, m, m);
92 delete [] work; 116
117 delete [] tau;
118 delete [] work;
119 }
93 } 120 }
94 121
95 /* 122 /*
96 ;;; Local Variables: *** 123 ;;; Local Variables: ***
97 ;;; mode: C++ *** 124 ;;; mode: C++ ***