Mercurial > hg > octave-nkf

comparison liboctave/CmplxQRP.cc @ 1922:8a57554f3142

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-11 02:53:20 by jwe]
author jwe
date Sun, 11 Feb 1996 02:53:20 +0000
parents 1281a23a34dd
children 20353fa5f83d
comparison
equal deleted inserted replaced
1921:ce0db0ca0729 1922:8a57554f3142
51 51
52 ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type) 52 ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type)
53 { 53 {
54 assert (qr_type != QR::raw); 54 assert (qr_type != QR::raw);
55 55
56 tau = 0;
57 work = 0;
58 tmp_data = 0;
59 jpvt = 0;
60
56 int m = a.rows (); 61 int m = a.rows ();
57 int n = a.cols (); 62 int n = a.cols ();
58 63
59 if (m == 0 || n == 0) 64 if (m == 0 || n == 0)
60 { 65 {
61 (*current_liboctave_error_handler) 66 (*current_liboctave_error_handler)
62 ("ComplexQR must have non-empty matrix"); 67 ("ComplexQR must have non-empty matrix");
63 return; 68 return;
64 } 69 }
65 70
66 Complex *tmp_data;
67 int min_mn = m < n ? m : n; 71 int min_mn = m < n ? m : n;
68 Complex *tau = new Complex[min_mn]; 72 tau = new Complex[min_mn];
73
69 int lwork = n; 74 int lwork = n;
70 Complex *work = new Complex[lwork]; 75 work = new Complex[lwork];
76
71 int info = 0; 77 int info = 0;
72 78
73 if (m > n) 79 if (m > n)
74 { 80 {
75 tmp_data = new Complex [m*m]; 81 tmp_data = new Complex [m*m];
76 copy (tmp_data, a.data (), a.length ()); 82 copy (tmp_data, a.data (), a.length ());
77 } 83 }
78 else 84 else
79 tmp_data = dup (a.data (), a.length ()); 85 tmp_data = dup (a.data (), a.length ());
80 86
81
82 work = new Complex[n]; 87 work = new Complex[n];
83 double *rwork = new double[2*n]; 88 rwork = new double[2*n];
84 int *jpvt = new int[n]; 89 jpvt = new int[n];
85 90
86 // Clear Pivot vector (code to enforce a certain permutation would 91 // Clear Pivot vector (code to enforce a certain permutation would
87 // go here...) 92 // go here...)
88 93
89 for (int i = 0; i < n; i++) 94 for (int i = 0; i < n; i++)
90 jpvt[i] = 0; 95 jpvt[i] = 0;
91 96
92 F77_FCN (zgeqpf, ZGEQPF) (m, n, tmp_data, m, jpvt, tau, work, rwork, 97 F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, jpvt, tau, work, rwork,
93 info); 98 info));
94 99
95 // Form Permutation matrix (if economy is requested, return the 100 delete [] work;
96 // indices only!) 101 work = 0;
97 102
98 if (qr_type == QR::economy && m > n) 103 delete [] rwork;
99 { 104 rwork = 0;
100 p.resize (1, n, 0.0); 105
101 for (int j = 0; j < n; j++) 106 if (f77_exception_encountered)
102 p.elem (0, j) = jpvt[j]; 107 (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf");
103 }
104 else 108 else
105 { 109 {
106 p.resize (n, n, 0.0); 110 // Form Permutation matrix (if economy is requested, return the
111 // indices only!)
112
113 if (qr_type == QR::economy && m > n)
114 {
115 p.resize (1, n, 0.0);
116 for (int j = 0; j < n; j++)
117 p.elem (0, j) = jpvt[j];
118 }
119 else
120 {
121 p.resize (n, n, 0.0);
122 for (int j = 0; j < n; j++)
123 p.elem (jpvt[j]-1, j) = 1.0;
124 }
125
126 delete [] jpvt;
127 jpvt = 0;
128
129 volatile int n2;
130
131 if (qr_type == QR::economy && m > n)
132 {
133 n2 = n;
134 r.resize (n, n, 0.0);
135 }
136 else
137 {
138 n2 = m;
139 r.resize (m, n, 0.0);
140 }
141
107 for (int j = 0; j < n; j++) 142 for (int j = 0; j < n; j++)
108 p.elem (jpvt[j]-1, j) = 1.0; 143 {
144 int limit = j < min_mn-1 ? j : min_mn-1;
145 for (int i = 0; i <= limit; i++)
146 r.elem (i, j) = tmp_data[m*j+i];
147 }
148
149 lwork = 32*m;
150 work = new Complex[lwork];
151
152 F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, tau, work,
153 lwork, info));
154
155 if (f77_exception_encountered)
156 (*current_liboctave_error_handler) ("unrecoverable error in zungqr");
157 else
158 {
159 q = ComplexMatrix (tmp_data, m, m);
160 tmp_data = 0;
161 q.resize (m, n2);
162 }
109 } 163 }
110 164
165 delete [] tau;
166 tau = 0;
167
168 delete [] work;
169 work = 0;
170
171 delete [] tmp_data;
172 tmp_data = 0;
173
111 delete [] jpvt; 174 delete [] jpvt;
112 delete [] work; 175 jpvt = 0;
113
114 int n2;
115 if (qr_type == QR::economy && m > n)
116 {
117 n2 = n;
118 r.resize (n, n, 0.0);
119 }
120 else
121 {
122 n2 = m;
123 r.resize (m, n, 0.0);
124 }
125
126 for (int j = 0; j < n; j++)
127 {
128 int limit = j < min_mn-1 ? j : min_mn-1;
129 for (int i = 0; i <= limit; i++)
130 r.elem (i, j) = tmp_data[m*j+i];
131 }
132
133 lwork = 32*m;
134 work = new Complex[lwork];
135
136 F77_FCN (zungqr, ZUNGQR) (m, m, min_mn, tmp_data, m, tau, work,
137 lwork, info);
138
139 q = ComplexMatrix (tmp_data, m, m);
140 q.resize (m, n2);
141
142 delete [] tau;
143 delete [] work;
144 } 176 }
145 177
146 /* 178 /*
147 ;;; Local Variables: *** 179 ;;; Local Variables: ***
148 ;;; mode: C++ *** 180 ;;; mode: C++ ***