Mercurial > hg > octave-nkf

comparison src/xdiv.cc @ 1800:024c75af53f1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-01-29 05:55:45 by jwe]
author jwe
date Mon, 29 Jan 1996 05:55:45 +0000
parents a02f140ed897
children e62277bf5fe0
comparison
equal deleted inserted replaced
1799:1881b2a7d1e2 1800:024c75af53f1
30 #include "CMatrix.h" 30 #include "CMatrix.h"
31 #include "dMatrix.h" 31 #include "dMatrix.h"
32 #include "oct-cmplx.h" 32 #include "oct-cmplx.h"
33 33
34 #include "error.h" 34 #include "error.h"
35 #include "pt-const.h"
36 #include "xdiv.h" 35 #include "xdiv.h"
37 36
38 static inline int 37 static inline int
39 result_ok (int info, double rcond, int warn = 1) 38 result_ok (int info, double rcond, int warn = 1)
40 { 39 {
85 // +---+----+ 84 // +---+----+
86 // complex_matrix | 2 | 4 | 85 // complex_matrix | 2 | 4 |
87 // +---+----+ 86 // +---+----+
88 87
89 // -*- 1 -*- 88 // -*- 1 -*-
90 tree_constant 89 Matrix
91 xdiv (const Matrix& a, const Matrix& b) 90 xdiv (const Matrix& a, const Matrix& b)
92 { 91 {
93 if (! mx_div_conform (b.columns (), a.columns ())) 92 if (! mx_div_conform (b.columns (), a.columns ()))
94 return tree_constant (); 93 return Matrix ();
95 94
96 Matrix atmp = a.transpose (); 95 Matrix atmp = a.transpose ();
97 Matrix btmp = b.transpose (); 96 Matrix btmp = b.transpose ();
98 97
99 int info; 98 int info;
100 if (btmp.rows () == btmp.columns ()) 99 if (btmp.rows () == btmp.columns ())
101 { 100 {
102 double rcond = 0.0; 101 double rcond = 0.0;
103 Matrix result = btmp.solve (atmp, info, rcond); 102 Matrix result = btmp.solve (atmp, info, rcond);
104 if (result_ok (info, rcond)) 103 if (result_ok (info, rcond))
105 return tree_constant (result.transpose ()); 104 return Matrix (result.transpose ());
106 } 105 }
107 106
108 int rank; 107 int rank;
109 Matrix result = btmp.lssolve (atmp, info, rank); 108 Matrix result = btmp.lssolve (atmp, info, rank);
110 109
111 return tree_constant (result.transpose ()); 110 return result.transpose ();
112 } 111 }
113 112
114 // -*- 2 -*- 113 // -*- 2 -*-
115 tree_constant 114 ComplexMatrix
116 xdiv (const Matrix& a, const ComplexMatrix& b) 115 xdiv (const Matrix& a, const ComplexMatrix& b)
117 { 116 {
118 if (! mx_div_conform (b.columns (), a.columns ())) 117 if (! mx_div_conform (b.columns (), a.columns ()))
119 return tree_constant (); 118 return ComplexMatrix ();
120 119
121 Matrix atmp = a.transpose (); 120 Matrix atmp = a.transpose ();
122 ComplexMatrix btmp = b.hermitian (); 121 ComplexMatrix btmp = b.hermitian ();
123 122
124 int info; 123 int info;
125 if (btmp.rows () == btmp.columns ()) 124 if (btmp.rows () == btmp.columns ())
126 { 125 {
127 double rcond = 0.0; 126 double rcond = 0.0;
128 ComplexMatrix result = btmp.solve (atmp, info, rcond); 127 ComplexMatrix result = btmp.solve (atmp, info, rcond);
129 if (result_ok (info, rcond)) 128 if (result_ok (info, rcond))
130 return tree_constant (result.hermitian ()); 129 return result.hermitian ();
131 } 130 }
132 131
133 int rank; 132 int rank;
134 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 133 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
135 134
136 return tree_constant (result.hermitian ()); 135 return result.hermitian ();
137 } 136 }
138 137
139 // -*- 3 -*- 138 // -*- 3 -*-
140 tree_constant 139 ComplexMatrix
141 xdiv (const ComplexMatrix& a, const Matrix& b) 140 xdiv (const ComplexMatrix& a, const Matrix& b)
142 { 141 {
143 if (! mx_div_conform (b.columns (), a.columns ())) 142 if (! mx_div_conform (b.columns (), a.columns ()))
144 return tree_constant (); 143 return ComplexMatrix ();
145 144
146 ComplexMatrix atmp = a.hermitian (); 145 ComplexMatrix atmp = a.hermitian ();
147 Matrix btmp = b.transpose (); 146 Matrix btmp = b.transpose ();
148 147
149 int info; 148 int info;
150 if (btmp.rows () == btmp.columns ()) 149 if (btmp.rows () == btmp.columns ())
151 { 150 {
152 double rcond = 0.0; 151 double rcond = 0.0;
153 ComplexMatrix result = btmp.solve (atmp, info, rcond); 152 ComplexMatrix result = btmp.solve (atmp, info, rcond);
154 if (result_ok (info, rcond)) 153 if (result_ok (info, rcond))
155 return tree_constant (result.hermitian ()); 154 return result.hermitian ();
156 } 155 }
157 156
158 int rank; 157 int rank;
159 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 158 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
160 159
161 return tree_constant (result.hermitian ()); 160 return result.hermitian ();
162 } 161 }
163 162
164 // -*- 4 -*- 163 // -*- 4 -*-
165 tree_constant 164 ComplexMatrix
166 xdiv (const ComplexMatrix& a, const ComplexMatrix& b) 165 xdiv (const ComplexMatrix& a, const ComplexMatrix& b)
167 { 166 {
168 if (! mx_div_conform (b.columns (), a.columns ())) 167 if (! mx_div_conform (b.columns (), a.columns ()))
169 return tree_constant (); 168 return ComplexMatrix ();
170 169
171 ComplexMatrix atmp = a.hermitian (); 170 ComplexMatrix atmp = a.hermitian ();
172 ComplexMatrix btmp = b.hermitian (); 171 ComplexMatrix btmp = b.hermitian ();
173 172
174 int info; 173 int info;
175 if (btmp.rows () == btmp.columns ()) 174 if (btmp.rows () == btmp.columns ())
176 { 175 {
177 double rcond = 0.0; 176 double rcond = 0.0;
178 ComplexMatrix result = btmp.solve (atmp, info, rcond); 177 ComplexMatrix result = btmp.solve (atmp, info, rcond);
179 if (result_ok (info, rcond)) 178 if (result_ok (info, rcond))
180 return tree_constant (result.hermitian ()); 179 return result.hermitian ();
181 } 180 }
182 181
183 int rank; 182 int rank;
184 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 183 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
185 184
186 return tree_constant (result.hermitian ()); 185 return result.hermitian ();
187 } 186 }
188 187
189 // Funny element by element division operations. 188 // Funny element by element division operations.
190 // 189 //
191 // op2 \ op1: s cs 190 // op2 \ op1: s cs
193 // matrix | 1 | 3 | 192 // matrix | 1 | 3 |
194 // +---+----+ 193 // +---+----+
195 // complex_matrix | 2 | 4 | 194 // complex_matrix | 2 | 4 |
196 // +---+----+ 195 // +---+----+
197 196
198 tree_constant 197 Matrix
199 x_el_div (double a, const Matrix& b) 198 x_el_div (double a, const Matrix& b)
200 { 199 {
201 int nr = b.rows (); 200 int nr = b.rows ();
202 int nc = b.columns (); 201 int nc = b.columns ();
203 202
205 204
206 for (int j = 0; j < nc; j++) 205 for (int j = 0; j < nc; j++)
207 for (int i = 0; i < nr; i++) 206 for (int i = 0; i < nr; i++)
208 result.elem (i, j) = a / b.elem (i, j); 207 result.elem (i, j) = a / b.elem (i, j);
209 208
210 return tree_constant (result); 209 return result;
211 } 210 }
212 211
213 tree_constant 212 ComplexMatrix
214 x_el_div (double a, const ComplexMatrix& b) 213 x_el_div (double a, const ComplexMatrix& b)
215 { 214 {
216 int nr = b.rows (); 215 int nr = b.rows ();
217 int nc = b.columns (); 216 int nc = b.columns ();
218 217
220 219
221 for (int j = 0; j < nc; j++) 220 for (int j = 0; j < nc; j++)
222 for (int i = 0; i < nr; i++) 221 for (int i = 0; i < nr; i++)
223 result.elem (i, j) = a / b.elem (i, j); 222 result.elem (i, j) = a / b.elem (i, j);
224 223
225 return tree_constant (result); 224 return result;
226 } 225 }
227 226
228 tree_constant 227 ComplexMatrix
229 x_el_div (const Complex a, const Matrix& b) 228 x_el_div (const Complex a, const Matrix& b)
230 { 229 {
231 int nr = b.rows (); 230 int nr = b.rows ();
232 int nc = b.columns (); 231 int nc = b.columns ();
233 232
235 234
236 for (int j = 0; j < nc; j++) 235 for (int j = 0; j < nc; j++)
237 for (int i = 0; i < nr; i++) 236 for (int i = 0; i < nr; i++)
238 result.elem (i, j) = a / b.elem (i, j); 237 result.elem (i, j) = a / b.elem (i, j);
239 238
240 return tree_constant (result); 239 return result;
241 } 240 }
242 241
243 tree_constant 242 ComplexMatrix
244 x_el_div (const Complex a, const ComplexMatrix& b) 243 x_el_div (const Complex a, const ComplexMatrix& b)
245 { 244 {
246 int nr = b.rows (); 245 int nr = b.rows ();
247 int nc = b.columns (); 246 int nc = b.columns ();
248 247
250 249
251 for (int j = 0; j < nc; j++) 250 for (int j = 0; j < nc; j++)
252 for (int i = 0; i < nr; i++) 251 for (int i = 0; i < nr; i++)
253 result.elem (i, j) = a / b.elem (i, j); 252 result.elem (i, j) = a / b.elem (i, j);
254 253
255 return tree_constant (result); 254 return result;
256 } 255 }
257 256
258 // Left division functions. 257 // Left division functions.
259 // 258 //
260 // op2 \ op1: m cm 259 // op2 \ op1: m cm
263 // +---+----+ 262 // +---+----+
264 // complex_matrix | 2 | 4 | 263 // complex_matrix | 2 | 4 |
265 // +---+----+ 264 // +---+----+
266 265
267 // -*- 1 -*- 266 // -*- 1 -*-
268 tree_constant 267 Matrix
269 xleftdiv (const Matrix& a, const Matrix& b) 268 xleftdiv (const Matrix& a, const Matrix& b)
270 { 269 {
271 if (! mx_leftdiv_conform (a.rows (), b.rows ())) 270 if (! mx_leftdiv_conform (a.rows (), b.rows ()))
272 return tree_constant (); 271 return Matrix ();
273 272
274 int info; 273 int info;
275 if (a.rows () == a.columns ()) 274 if (a.rows () == a.columns ())
276 { 275 {
277 double rcond = 0.0; 276 double rcond = 0.0;
278 Matrix result = a.solve (b, info, rcond); 277 Matrix result = a.solve (b, info, rcond);
279 if (result_ok (info, rcond)) 278 if (result_ok (info, rcond))
280 return tree_constant (result); 279 return result;
281 } 280 }
282 281
283 int rank; 282 int rank;
284 Matrix result = a.lssolve (b, info, rank); 283 return a.lssolve (b, info, rank);
285
286 return tree_constant (result);
287 } 284 }
288 285
289 // -*- 2 -*- 286 // -*- 2 -*-
290 tree_constant 287 ComplexMatrix
291 xleftdiv (const Matrix& a, const ComplexMatrix& b) 288 xleftdiv (const Matrix& a, const ComplexMatrix& b)
292 { 289 {
293 if (! mx_leftdiv_conform (a.rows (), b.rows ())) 290 if (! mx_leftdiv_conform (a.rows (), b.rows ()))
294 return tree_constant (); 291 return ComplexMatrix ();
295 292
296 int info; 293 int info;
297 if (a.rows () == a.columns ()) 294 if (a.rows () == a.columns ())
298 { 295 {
299 double rcond = 0.0; 296 double rcond = 0.0;
300 ComplexMatrix result = a.solve (b, info, rcond); 297 ComplexMatrix result = a.solve (b, info, rcond);
301 if (result_ok (info, rcond)) 298 if (result_ok (info, rcond))
302 return tree_constant (result); 299 return result;
303 } 300 }
304 301
305 int rank; 302 int rank;
306 ComplexMatrix result = a.lssolve (b, info, rank); 303 return a.lssolve (b, info, rank);
307
308 return tree_constant (result);
309 } 304 }
310 305
311 // -*- 3 -*- 306 // -*- 3 -*-
312 tree_constant 307 ComplexMatrix
313 xleftdiv (const ComplexMatrix& a, const Matrix& b) 308 xleftdiv (const ComplexMatrix& a, const Matrix& b)
314 { 309 {
315 if (! mx_leftdiv_conform (a.rows (), b.rows ())) 310 if (! mx_leftdiv_conform (a.rows (), b.rows ()))
316 return tree_constant (); 311 return ComplexMatrix ();
317 312
318 int info; 313 int info;
319 if (a.rows () == a.columns ()) 314 if (a.rows () == a.columns ())
320 { 315 {
321 double rcond = 0.0; 316 double rcond = 0.0;
322 ComplexMatrix result = a.solve (b, info, rcond); 317 ComplexMatrix result = a.solve (b, info, rcond);
323 if (result_ok (info, rcond)) 318 if (result_ok (info, rcond))
324 return tree_constant (result); 319 return result;
325 } 320 }
326 321
327 int rank; 322 int rank;
328 ComplexMatrix result = a.lssolve (b, info, rank); 323 return a.lssolve (b, info, rank);
329
330 return tree_constant (result);
331 } 324 }
332 325
333 // -*- 4 -*- 326 // -*- 4 -*-
334 tree_constant 327 ComplexMatrix
335 xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b) 328 xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b)
336 { 329 {
337 if (! mx_leftdiv_conform (a.rows (), b.rows ())) 330 if (! mx_leftdiv_conform (a.rows (), b.rows ()))
338 return tree_constant (); 331 return ComplexMatrix ();
339 332
340 int info; 333 int info;
341 if (a.rows () == a.columns ()) 334 if (a.rows () == a.columns ())
342 { 335 {
343 double rcond = 0.0; 336 double rcond = 0.0;
344 ComplexMatrix result = a.solve (b, info, rcond); 337 ComplexMatrix result = a.solve (b, info, rcond);
345 if (result_ok (info, rcond)) 338 if (result_ok (info, rcond))
346 return tree_constant (result); 339 return result;
347 } 340 }
348 341
349 int rank; 342 int rank;
350 ComplexMatrix result = a.lssolve (b, info, rank); 343 return a.lssolve (b, info, rank);
351
352 return tree_constant (result);
353 } 344 }
354 345
355 /* 346 /*
356 ;;; Local Variables: *** 347 ;;; Local Variables: ***
357 ;;; mode: C++ *** 348 ;;; mode: C++ ***