comparison src/sparse-xdiv.cc @ 5275:23b37da9fd5b

[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents 57077d0ddc8e
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
34 #include "CSparse.h" 34 #include "CSparse.h"
35 #include "oct-spparms.h" 35 #include "oct-spparms.h"
36 #include "sparse-xdiv.h" 36 #include "sparse-xdiv.h"
37 37
38 static inline bool 38 static inline bool
39 result_ok (int info) 39 result_ok (octave_idx_type info)
40 { 40 {
41 #ifdef HAVE_LSSOLVE 41 #ifdef HAVE_LSSOLVE
42 return (info != -2 && info != -1); 42 return (info != -2 && info != -1);
43 #else 43 #else
44 // If the matrix is singular, who cares as we don't have QR based solver yet 44 // If the matrix is singular, who cares as we don't have QR based solver yet
55 55
56 template <class T1, class T2> 56 template <class T1, class T2>
57 bool 57 bool
58 mx_leftdiv_conform (const T1& a, const T2& b) 58 mx_leftdiv_conform (const T1& a, const T2& b)
59 { 59 {
60 int a_nr = a.rows (); 60 octave_idx_type a_nr = a.rows ();
61 int b_nr = b.rows (); 61 octave_idx_type b_nr = b.rows ();
62 62
63 if (a_nr != b_nr) 63 if (a_nr != b_nr)
64 { 64 {
65 int a_nc = a.cols (); 65 octave_idx_type a_nc = a.cols ();
66 int b_nc = b.cols (); 66 octave_idx_type b_nc = b.cols ();
67 67
68 gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); 68 gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
69 return false; 69 return false;
70 } 70 }
71 71
86 86
87 template <class T1, class T2> 87 template <class T1, class T2>
88 bool 88 bool
89 mx_div_conform (const T1& a, const T2& b) 89 mx_div_conform (const T1& a, const T2& b)
90 { 90 {
91 int a_nc = a.cols (); 91 octave_idx_type a_nc = a.cols ();
92 int b_nc = b.cols (); 92 octave_idx_type b_nc = b.cols ();
93 93
94 if (a_nc != b_nc) 94 if (a_nc != b_nc)
95 { 95 {
96 int a_nr = a.rows (); 96 octave_idx_type a_nr = a.rows ();
97 int b_nr = b.rows (); 97 octave_idx_type b_nr = b.rows ();
98 98
99 gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); 99 gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
100 return false; 100 return false;
101 } 101 }
102 102
132 return Matrix (); 132 return Matrix ();
133 133
134 Matrix atmp = a.transpose (); 134 Matrix atmp = a.transpose ();
135 SparseMatrix btmp = b.transpose (); 135 SparseMatrix btmp = b.transpose ();
136 136
137 int info; 137 octave_idx_type info;
138 if (btmp.rows () == btmp.columns ()) 138 if (btmp.rows () == btmp.columns ())
139 { 139 {
140 double rcond = 0.0; 140 double rcond = 0.0;
141 141
142 Matrix result = btmp.solve (atmp, info, rcond, 142 Matrix result = btmp.solve (atmp, info, rcond,
144 144
145 if (result_ok (info)) 145 if (result_ok (info))
146 return Matrix (result.transpose ()); 146 return Matrix (result.transpose ());
147 } 147 }
148 148
149 int rank; 149 octave_idx_type rank;
150 Matrix result = btmp.lssolve (atmp, info, rank); 150 Matrix result = btmp.lssolve (atmp, info, rank);
151 151
152 return result.transpose (); 152 return result.transpose ();
153 } 153 }
154 154
160 return ComplexMatrix (); 160 return ComplexMatrix ();
161 161
162 Matrix atmp = a.transpose (); 162 Matrix atmp = a.transpose ();
163 SparseComplexMatrix btmp = b.hermitian (); 163 SparseComplexMatrix btmp = b.hermitian ();
164 164
165 int info; 165 octave_idx_type info;
166 if (btmp.rows () == btmp.columns ()) 166 if (btmp.rows () == btmp.columns ())
167 { 167 {
168 double rcond = 0.0; 168 double rcond = 0.0;
169 169
170 ComplexMatrix result 170 ComplexMatrix result
172 172
173 if (result_ok (info)) 173 if (result_ok (info))
174 return result.hermitian (); 174 return result.hermitian ();
175 } 175 }
176 176
177 int rank; 177 octave_idx_type rank;
178 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 178 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
179 179
180 return result.hermitian (); 180 return result.hermitian ();
181 } 181 }
182 182
188 return ComplexMatrix (); 188 return ComplexMatrix ();
189 189
190 ComplexMatrix atmp = a.hermitian (); 190 ComplexMatrix atmp = a.hermitian ();
191 SparseMatrix btmp = b.transpose (); 191 SparseMatrix btmp = b.transpose ();
192 192
193 int info; 193 octave_idx_type info;
194 if (btmp.rows () == btmp.columns ()) 194 if (btmp.rows () == btmp.columns ())
195 { 195 {
196 double rcond = 0.0; 196 double rcond = 0.0;
197 197
198 ComplexMatrix result 198 ComplexMatrix result
200 200
201 if (result_ok (info)) 201 if (result_ok (info))
202 return result.hermitian (); 202 return result.hermitian ();
203 } 203 }
204 204
205 int rank; 205 octave_idx_type rank;
206 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 206 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
207 207
208 return result.hermitian (); 208 return result.hermitian ();
209 } 209 }
210 210
216 return ComplexMatrix (); 216 return ComplexMatrix ();
217 217
218 ComplexMatrix atmp = a.hermitian (); 218 ComplexMatrix atmp = a.hermitian ();
219 SparseComplexMatrix btmp = b.hermitian (); 219 SparseComplexMatrix btmp = b.hermitian ();
220 220
221 int info; 221 octave_idx_type info;
222 if (btmp.rows () == btmp.columns ()) 222 if (btmp.rows () == btmp.columns ())
223 { 223 {
224 double rcond = 0.0; 224 double rcond = 0.0;
225 225
226 ComplexMatrix result 226 ComplexMatrix result
228 228
229 if (result_ok (info)) 229 if (result_ok (info))
230 return result.hermitian (); 230 return result.hermitian ();
231 } 231 }
232 232
233 int rank; 233 octave_idx_type rank;
234 ComplexMatrix result = btmp.lssolve (atmp, info, rank); 234 ComplexMatrix result = btmp.lssolve (atmp, info, rank);
235 235
236 return result.hermitian (); 236 return result.hermitian ();
237 } 237 }
238 238
244 return SparseMatrix (); 244 return SparseMatrix ();
245 245
246 SparseMatrix atmp = a.transpose (); 246 SparseMatrix atmp = a.transpose ();
247 SparseMatrix btmp = b.transpose (); 247 SparseMatrix btmp = b.transpose ();
248 248
249 int info; 249 octave_idx_type info;
250 if (btmp.rows () == btmp.columns ()) 250 if (btmp.rows () == btmp.columns ())
251 { 251 {
252 double rcond = 0.0; 252 double rcond = 0.0;
253 253
254 SparseMatrix result = btmp.solve (atmp, info, rcond, 254 SparseMatrix result = btmp.solve (atmp, info, rcond,
256 256
257 if (result_ok (info)) 257 if (result_ok (info))
258 return SparseMatrix (result.transpose ()); 258 return SparseMatrix (result.transpose ());
259 } 259 }
260 260
261 int rank; 261 octave_idx_type rank;
262 SparseMatrix result = btmp.lssolve (atmp, info, rank); 262 SparseMatrix result = btmp.lssolve (atmp, info, rank);
263 263
264 return result.transpose (); 264 return result.transpose ();
265 } 265 }
266 266
272 return SparseComplexMatrix (); 272 return SparseComplexMatrix ();
273 273
274 SparseMatrix atmp = a.transpose (); 274 SparseMatrix atmp = a.transpose ();
275 SparseComplexMatrix btmp = b.hermitian (); 275 SparseComplexMatrix btmp = b.hermitian ();
276 276
277 int info; 277 octave_idx_type info;
278 if (btmp.rows () == btmp.columns ()) 278 if (btmp.rows () == btmp.columns ())
279 { 279 {
280 double rcond = 0.0; 280 double rcond = 0.0;
281 281
282 SparseComplexMatrix result 282 SparseComplexMatrix result
284 284
285 if (result_ok (info)) 285 if (result_ok (info))
286 return result.hermitian (); 286 return result.hermitian ();
287 } 287 }
288 288
289 int rank; 289 octave_idx_type rank;
290 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); 290 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
291 291
292 return result.hermitian (); 292 return result.hermitian ();
293 } 293 }
294 294
300 return SparseComplexMatrix (); 300 return SparseComplexMatrix ();
301 301
302 SparseComplexMatrix atmp = a.hermitian (); 302 SparseComplexMatrix atmp = a.hermitian ();
303 SparseMatrix btmp = b.transpose (); 303 SparseMatrix btmp = b.transpose ();
304 304
305 int info; 305 octave_idx_type info;
306 if (btmp.rows () == btmp.columns ()) 306 if (btmp.rows () == btmp.columns ())
307 { 307 {
308 double rcond = 0.0; 308 double rcond = 0.0;
309 309
310 SparseComplexMatrix result 310 SparseComplexMatrix result
312 312
313 if (result_ok (info)) 313 if (result_ok (info))
314 return result.hermitian (); 314 return result.hermitian ();
315 } 315 }
316 316
317 int rank; 317 octave_idx_type rank;
318 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); 318 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
319 319
320 return result.hermitian (); 320 return result.hermitian ();
321 } 321 }
322 322
328 return SparseComplexMatrix (); 328 return SparseComplexMatrix ();
329 329
330 SparseComplexMatrix atmp = a.hermitian (); 330 SparseComplexMatrix atmp = a.hermitian ();
331 SparseComplexMatrix btmp = b.hermitian (); 331 SparseComplexMatrix btmp = b.hermitian ();
332 332
333 int info; 333 octave_idx_type info;
334 if (btmp.rows () == btmp.columns ()) 334 if (btmp.rows () == btmp.columns ())
335 { 335 {
336 double rcond = 0.0; 336 double rcond = 0.0;
337 337
338 SparseComplexMatrix result 338 SparseComplexMatrix result
340 340
341 if (result_ok (info)) 341 if (result_ok (info))
342 return result.hermitian (); 342 return result.hermitian ();
343 } 343 }
344 344
345 int rank; 345 octave_idx_type rank;
346 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); 346 SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
347 347
348 return result.hermitian (); 348 return result.hermitian ();
349 } 349 }
350 350
358 // +---+----+ 358 // +---+----+
359 359
360 Matrix 360 Matrix
361 x_el_div (double a, const SparseMatrix& b) 361 x_el_div (double a, const SparseMatrix& b)
362 { 362 {
363 int nr = b.rows (); 363 octave_idx_type nr = b.rows ();
364 int nc = b.cols (); 364 octave_idx_type nc = b.cols ();
365 365
366 Matrix result; 366 Matrix result;
367 if (a == 0.) 367 if (a == 0.)
368 result = Matrix (nr, nc, octave_NaN); 368 result = Matrix (nr, nc, octave_NaN);
369 else if (a > 0.) 369 else if (a > 0.)
370 result = Matrix (nr, nc, octave_Inf); 370 result = Matrix (nr, nc, octave_Inf);
371 else 371 else
372 result = Matrix (nr, nc, -octave_Inf); 372 result = Matrix (nr, nc, -octave_Inf);
373 373
374 374
375 for (int j = 0; j < nc; j++) 375 for (octave_idx_type j = 0; j < nc; j++)
376 for (int i = b.cidx(j); i < b.cidx(j+1); i++) 376 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
377 { 377 {
378 OCTAVE_QUIT; 378 OCTAVE_QUIT;
379 result.elem (b.ridx(i), j) = a / b.data (i); 379 result.elem (b.ridx(i), j) = a / b.data (i);
380 } 380 }
381 381
383 } 383 }
384 384
385 ComplexMatrix 385 ComplexMatrix
386 x_el_div (double a, const SparseComplexMatrix& b) 386 x_el_div (double a, const SparseComplexMatrix& b)
387 { 387 {
388 int nr = b.rows (); 388 octave_idx_type nr = b.rows ();
389 int nc = b.cols (); 389 octave_idx_type nc = b.cols ();
390 390
391 ComplexMatrix result (nr, nc, Complex(octave_NaN, octave_NaN)); 391 ComplexMatrix result (nr, nc, Complex(octave_NaN, octave_NaN));
392 392
393 for (int j = 0; j < nc; j++) 393 for (octave_idx_type j = 0; j < nc; j++)
394 for (int i = b.cidx(j); i < b.cidx(j+1); i++) 394 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
395 { 395 {
396 OCTAVE_QUIT; 396 OCTAVE_QUIT;
397 result.elem (b.ridx(i), j) = a / b.data (i); 397 result.elem (b.ridx(i), j) = a / b.data (i);
398 } 398 }
399 399
401 } 401 }
402 402
403 ComplexMatrix 403 ComplexMatrix
404 x_el_div (const Complex a, const SparseMatrix& b) 404 x_el_div (const Complex a, const SparseMatrix& b)
405 { 405 {
406 int nr = b.rows (); 406 octave_idx_type nr = b.rows ();
407 int nc = b.cols (); 407 octave_idx_type nc = b.cols ();
408 408
409 ComplexMatrix result (nr, nc, (a / 0.0)); 409 ComplexMatrix result (nr, nc, (a / 0.0));
410 410
411 for (int j = 0; j < nc; j++) 411 for (octave_idx_type j = 0; j < nc; j++)
412 for (int i = b.cidx(j); i < b.cidx(j+1); i++) 412 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
413 { 413 {
414 OCTAVE_QUIT; 414 OCTAVE_QUIT;
415 result.elem (b.ridx(i), j) = a / b.data (i); 415 result.elem (b.ridx(i), j) = a / b.data (i);
416 } 416 }
417 417
419 } 419 }
420 420
421 ComplexMatrix 421 ComplexMatrix
422 x_el_div (const Complex a, const SparseComplexMatrix& b) 422 x_el_div (const Complex a, const SparseComplexMatrix& b)
423 { 423 {
424 int nr = b.rows (); 424 octave_idx_type nr = b.rows ();
425 int nc = b.cols (); 425 octave_idx_type nc = b.cols ();
426 426
427 ComplexMatrix result (nr, nc, (a / 0.0)); 427 ComplexMatrix result (nr, nc, (a / 0.0));
428 428
429 for (int j = 0; j < nc; j++) 429 for (octave_idx_type j = 0; j < nc; j++)
430 for (int i = b.cidx(j); i < b.cidx(j+1); i++) 430 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
431 { 431 {
432 OCTAVE_QUIT; 432 OCTAVE_QUIT;
433 result.elem (b.ridx(i), j) = a / b.data (i); 433 result.elem (b.ridx(i), j) = a / b.data (i);
434 } 434 }
435 435
454 xleftdiv (const SparseMatrix& a, const Matrix& b) 454 xleftdiv (const SparseMatrix& a, const Matrix& b)
455 { 455 {
456 if (! mx_leftdiv_conform (a, b)) 456 if (! mx_leftdiv_conform (a, b))
457 return Matrix (); 457 return Matrix ();
458 458
459 int info; 459 octave_idx_type info;
460 if (a.rows () == a.columns ()) 460 if (a.rows () == a.columns ())
461 { 461 {
462 double rcond = 0.0; 462 double rcond = 0.0;
463 463
464 Matrix result 464 Matrix result
466 466
467 if (result_ok (info)) 467 if (result_ok (info))
468 return result; 468 return result;
469 } 469 }
470 470
471 int rank; 471 octave_idx_type rank;
472 return a.lssolve (b, info, rank); 472 return a.lssolve (b, info, rank);
473 } 473 }
474 474
475 // -*- 2 -*- 475 // -*- 2 -*-
476 ComplexMatrix 476 ComplexMatrix
477 xleftdiv (const SparseMatrix& a, const ComplexMatrix& b) 477 xleftdiv (const SparseMatrix& a, const ComplexMatrix& b)
478 { 478 {
479 if (! mx_leftdiv_conform (a, b)) 479 if (! mx_leftdiv_conform (a, b))
480 return ComplexMatrix (); 480 return ComplexMatrix ();
481 481
482 int info; 482 octave_idx_type info;
483 if (a.rows () == a.columns ()) 483 if (a.rows () == a.columns ())
484 { 484 {
485 double rcond = 0.0; 485 double rcond = 0.0;
486 486
487 ComplexMatrix result 487 ComplexMatrix result
489 489
490 if (result_ok (info)) 490 if (result_ok (info))
491 return result; 491 return result;
492 } 492 }
493 493
494 int rank; 494 octave_idx_type rank;
495 return a.lssolve (b, info, rank); 495 return a.lssolve (b, info, rank);
496 } 496 }
497 497
498 // -*- 3 -*- 498 // -*- 3 -*-
499 SparseMatrix 499 SparseMatrix
500 xleftdiv (const SparseMatrix& a, const SparseMatrix& b) 500 xleftdiv (const SparseMatrix& a, const SparseMatrix& b)
501 { 501 {
502 if (! mx_leftdiv_conform (a, b)) 502 if (! mx_leftdiv_conform (a, b))
503 return SparseMatrix (); 503 return SparseMatrix ();
504 504
505 int info; 505 octave_idx_type info;
506 if (a.rows () == a.columns ()) 506 if (a.rows () == a.columns ())
507 { 507 {
508 double rcond = 0.0; 508 double rcond = 0.0;
509 509
510 SparseMatrix result 510 SparseMatrix result
512 512
513 if (result_ok (info)) 513 if (result_ok (info))
514 return result; 514 return result;
515 } 515 }
516 516
517 int rank; 517 octave_idx_type rank;
518 return a.lssolve (b, info, rank); 518 return a.lssolve (b, info, rank);
519 } 519 }
520 520
521 // -*- 4 -*- 521 // -*- 4 -*-
522 SparseComplexMatrix 522 SparseComplexMatrix
523 xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b) 523 xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b)
524 { 524 {
525 if (! mx_leftdiv_conform (a, b)) 525 if (! mx_leftdiv_conform (a, b))
526 return SparseComplexMatrix (); 526 return SparseComplexMatrix ();
527 527
528 int info; 528 octave_idx_type info;
529 if (a.rows () == a.columns ()) 529 if (a.rows () == a.columns ())
530 { 530 {
531 double rcond = 0.0; 531 double rcond = 0.0;
532 532
533 SparseComplexMatrix result 533 SparseComplexMatrix result
535 535
536 if (result_ok (info)) 536 if (result_ok (info))
537 return result; 537 return result;
538 } 538 }
539 539
540 int rank; 540 octave_idx_type rank;
541 return a.lssolve (b, info, rank); 541 return a.lssolve (b, info, rank);
542 } 542 }
543 543
544 // -*- 5 -*- 544 // -*- 5 -*-
545 ComplexMatrix 545 ComplexMatrix
546 xleftdiv (const SparseComplexMatrix& a, const Matrix& b) 546 xleftdiv (const SparseComplexMatrix& a, const Matrix& b)
547 { 547 {
548 if (! mx_leftdiv_conform (a, b)) 548 if (! mx_leftdiv_conform (a, b))
549 return ComplexMatrix (); 549 return ComplexMatrix ();
550 550
551 int info; 551 octave_idx_type info;
552 if (a.rows () == a.columns ()) 552 if (a.rows () == a.columns ())
553 { 553 {
554 double rcond = 0.0; 554 double rcond = 0.0;
555 555
556 ComplexMatrix result 556 ComplexMatrix result
558 558
559 if (result_ok (info)) 559 if (result_ok (info))
560 return result; 560 return result;
561 } 561 }
562 562
563 int rank; 563 octave_idx_type rank;
564 return a.lssolve (b, info, rank); 564 return a.lssolve (b, info, rank);
565 } 565 }
566 566
567 // -*- 6 -*- 567 // -*- 6 -*-
568 ComplexMatrix 568 ComplexMatrix
569 xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b) 569 xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b)
570 { 570 {
571 if (! mx_leftdiv_conform (a, b)) 571 if (! mx_leftdiv_conform (a, b))
572 return ComplexMatrix (); 572 return ComplexMatrix ();
573 573
574 int info; 574 octave_idx_type info;
575 if (a.rows () == a.columns ()) 575 if (a.rows () == a.columns ())
576 { 576 {
577 double rcond = 0.0; 577 double rcond = 0.0;
578 578
579 ComplexMatrix result 579 ComplexMatrix result
581 581
582 if (result_ok (info)) 582 if (result_ok (info))
583 return result; 583 return result;
584 } 584 }
585 585
586 int rank; 586 octave_idx_type rank;
587 return a.lssolve (b, info, rank); 587 return a.lssolve (b, info, rank);
588 } 588 }
589 589
590 // -*- 7 -*- 590 // -*- 7 -*-
591 SparseComplexMatrix 591 SparseComplexMatrix
592 xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b) 592 xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b)
593 { 593 {
594 if (! mx_leftdiv_conform (a, b)) 594 if (! mx_leftdiv_conform (a, b))
595 return SparseComplexMatrix (); 595 return SparseComplexMatrix ();
596 596
597 int info; 597 octave_idx_type info;
598 if (a.rows () == a.columns ()) 598 if (a.rows () == a.columns ())
599 { 599 {
600 double rcond = 0.0; 600 double rcond = 0.0;
601 601
602 SparseComplexMatrix result 602 SparseComplexMatrix result
604 604
605 if (result_ok (info)) 605 if (result_ok (info))
606 return result; 606 return result;
607 } 607 }
608 608
609 int rank; 609 octave_idx_type rank;
610 return a.lssolve (b, info, rank); 610 return a.lssolve (b, info, rank);
611 } 611 }
612 612
613 // -*- 8 -*- 613 // -*- 8 -*-
614 SparseComplexMatrix 614 SparseComplexMatrix
615 xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b) 615 xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
616 { 616 {
617 if (! mx_leftdiv_conform (a, b)) 617 if (! mx_leftdiv_conform (a, b))
618 return SparseComplexMatrix (); 618 return SparseComplexMatrix ();
619 619
620 int info; 620 octave_idx_type info;
621 if (a.rows () == a.columns ()) 621 if (a.rows () == a.columns ())
622 { 622 {
623 double rcond = 0.0; 623 double rcond = 0.0;
624 624
625 SparseComplexMatrix result 625 SparseComplexMatrix result
627 627
628 if (result_ok (info)) 628 if (result_ok (info))
629 return result; 629 return result;
630 } 630 }
631 631
632 int rank; 632 octave_idx_type rank;
633 return a.lssolve (b, info, rank); 633 return a.lssolve (b, info, rank);
634 } 634 }
635 635
636 /* 636 /*
637 ;;; Local Variables: *** 637 ;;; Local Variables: ***