Mercurial > octave-nkf

comparison liboctave/EIG.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents 57632dea2446
children 72c96de7a403
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
71 octave_idx_type& 71 octave_idx_type&
72 F77_CHAR_ARG_LEN_DECL 72 F77_CHAR_ARG_LEN_DECL
73 F77_CHAR_ARG_LEN_DECL); 73 F77_CHAR_ARG_LEN_DECL);
74 74
75 F77_RET_T 75 F77_RET_T
76 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, 76 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL,
77 const octave_idx_type&, double*, 77 const octave_idx_type&, double*,
78 const octave_idx_type&, octave_idx_type& 78 const octave_idx_type&, octave_idx_type&
79 F77_CHAR_ARG_LEN_DECL 79 F77_CHAR_ARG_LEN_DECL
80 F77_CHAR_ARG_LEN_DECL); 80 F77_CHAR_ARG_LEN_DECL);
81 81
82 F77_RET_T 82 F77_RET_T
83 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, 83 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL,
84 const octave_idx_type&, 84 const octave_idx_type&,
85 Complex*, const octave_idx_type&, 85 Complex*, const octave_idx_type&,
86 octave_idx_type& 86 octave_idx_type&
87 F77_CHAR_ARG_LEN_DECL 87 F77_CHAR_ARG_LEN_DECL
88 F77_CHAR_ARG_LEN_DECL); 88 F77_CHAR_ARG_LEN_DECL);
89 89
90 F77_RET_T 90 F77_RET_T
91 F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL, 91 F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL,
92 F77_CONST_CHAR_ARG_DECL, 92 F77_CONST_CHAR_ARG_DECL,
93 const octave_idx_type&, 93 const octave_idx_type&,
94 double*, const octave_idx_type&, 94 double*, const octave_idx_type&,
95 double*, const octave_idx_type&, 95 double*, const octave_idx_type&,
96 double*, double*, double *, double*, 96 double*, double*, double *, double*,
97 const octave_idx_type&, double*, 97 const octave_idx_type&, double*,
98 const octave_idx_type&, double*, 98 const octave_idx_type&, double*,
110 const octave_idx_type&, octave_idx_type& 110 const octave_idx_type&, octave_idx_type&
111 F77_CHAR_ARG_LEN_DECL 111 F77_CHAR_ARG_LEN_DECL
112 F77_CHAR_ARG_LEN_DECL); 112 F77_CHAR_ARG_LEN_DECL);
113 113
114 F77_RET_T 114 F77_RET_T
115 F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL, 115 F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL,
116 F77_CONST_CHAR_ARG_DECL, 116 F77_CONST_CHAR_ARG_DECL,
117 const octave_idx_type&, 117 const octave_idx_type&,
118 Complex*, const octave_idx_type&, 118 Complex*, const octave_idx_type&,
119 Complex*, const octave_idx_type&, 119 Complex*, const octave_idx_type&,
120 Complex*, Complex*, Complex*, 120 Complex*, Complex*, Complex*,
121 const octave_idx_type&, Complex*, 121 const octave_idx_type&, Complex*,
122 const octave_idx_type&, Complex*, 122 const octave_idx_type&, Complex*,
124 F77_CHAR_ARG_LEN_DECL 124 F77_CHAR_ARG_LEN_DECL
125 F77_CHAR_ARG_LEN_DECL); 125 F77_CHAR_ARG_LEN_DECL);
126 126
127 F77_RET_T 127 F77_RET_T
128 F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&, 128 F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&,
129 F77_CONST_CHAR_ARG_DECL, 129 F77_CONST_CHAR_ARG_DECL,
130 F77_CONST_CHAR_ARG_DECL, 130 F77_CONST_CHAR_ARG_DECL,
131 const octave_idx_type&, Complex*, 131 const octave_idx_type&, Complex*,
132 const octave_idx_type&, Complex*, 132 const octave_idx_type&, Complex*,
133 const octave_idx_type&, double*, Complex*, 133 const octave_idx_type&, double*, Complex*,
134 const octave_idx_type&, double*, octave_idx_type& 134 const octave_idx_type&, double*, octave_idx_type&
248 (*current_liboctave_error_handler) ("dgeev workspace query failed"); 248 (*current_liboctave_error_handler) ("dgeev workspace query failed");
249 249
250 return info; 250 return info;
251 } 251 }
252 252
253 octave_idx_type 253 octave_idx_type
254 EIG::symmetric_init (const Matrix& a, bool calc_ev) 254 EIG::symmetric_init (const Matrix& a, bool calc_ev)
255 { 255 {
256 octave_idx_type n = a.rows (); 256 octave_idx_type n = a.rows ();
257 257
258 if (n != a.cols ()) 258 if (n != a.cols ())
490 490
491 Matrix tmp = b; 491 Matrix tmp = b;
492 double *tmp_data = tmp.fortran_vec (); 492 double *tmp_data = tmp.fortran_vec ();
493 493
494 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), 494 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
495 n, tmp_data, n, 495 n, tmp_data, n,
496 info 496 info
497 F77_CHAR_ARG_LEN (1) 497 F77_CHAR_ARG_LEN (1)
498 F77_CHAR_ARG_LEN (1))); 498 F77_CHAR_ARG_LEN (1)));
499 499
500 if (a.is_symmetric () && b.is_symmetric () && info == 0) 500 if (a.is_symmetric () && b.is_symmetric () && info == 0)
525 double *dummy = 0; 525 double *dummy = 0;
526 octave_idx_type idummy = 1; 526 octave_idx_type idummy = 1;
527 527
528 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), 528 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
529 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 529 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
530 n, atmp_data, n, btmp_data, n, 530 n, atmp_data, n, btmp_data, n,
531 par, pai, pbeta, 531 par, pai, pbeta,
532 dummy, idummy, pvr, n, 532 dummy, idummy, pvr, n,
533 &dummy_work, lwork, info 533 &dummy_work, lwork, info
534 F77_CHAR_ARG_LEN (1) 534 F77_CHAR_ARG_LEN (1)
535 F77_CHAR_ARG_LEN (1))); 535 F77_CHAR_ARG_LEN (1)));
540 Array<double> work (dim_vector (lwork, 1)); 540 Array<double> work (dim_vector (lwork, 1));
541 double *pwork = work.fortran_vec (); 541 double *pwork = work.fortran_vec ();
542 542
543 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), 543 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
544 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 544 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
545 n, atmp_data, n, btmp_data, n, 545 n, atmp_data, n, btmp_data, n,
546 par, pai, pbeta, 546 par, pai, pbeta,
547 dummy, idummy, pvr, n, 547 dummy, idummy, pvr, n,
548 pwork, lwork, info 548 pwork, lwork, info
549 F77_CHAR_ARG_LEN (1) 549 F77_CHAR_ARG_LEN (1)
550 F77_CHAR_ARG_LEN (1))); 550 F77_CHAR_ARG_LEN (1)));
579 { 579 {
580 (*current_liboctave_error_handler) ("EIG: internal error"); 580 (*current_liboctave_error_handler) ("EIG: internal error");
581 return -1; 581 return -1;
582 } 582 }
583 583
584 lambda.elem(j) = Complex (ar.elem(j) / beta.elem (j), 584 lambda.elem(j) = Complex (ar.elem(j) / beta.elem (j),
585 ai.elem(j) / beta.elem (j)); 585 ai.elem(j) / beta.elem (j));
586 lambda.elem(j+1) = Complex (ar.elem(j+1) / beta.elem (j+1), 586 lambda.elem(j+1) = Complex (ar.elem(j+1) / beta.elem (j+1),
587 ai.elem(j+1) / beta.elem (j+1)); 587 ai.elem(j+1) / beta.elem (j+1));
588 588
589 for (octave_idx_type i = 0; i < nvr; i++) 589 for (octave_idx_type i = 0; i < nvr; i++)
590 { 590 {
591 double real_part = vr.elem (i, j); 591 double real_part = vr.elem (i, j);
601 (*current_liboctave_error_handler) ("dggev workspace query failed"); 601 (*current_liboctave_error_handler) ("dggev workspace query failed");
602 602
603 return info; 603 return info;
604 } 604 }
605 605
606 octave_idx_type 606 octave_idx_type
607 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev) 607 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev)
608 { 608 {
609 octave_idx_type n = a.rows (); 609 octave_idx_type n = a.rows ();
610 octave_idx_type nb = b.rows (); 610 octave_idx_type nb = b.rows ();
611 611
635 octave_idx_type lwork = -1; 635 octave_idx_type lwork = -1;
636 double dummy_work; 636 double dummy_work;
637 637
638 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 638 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
639 F77_CONST_CHAR_ARG2 ("U", 1), 639 F77_CONST_CHAR_ARG2 ("U", 1),
640 n, atmp_data, n, 640 n, atmp_data, n,
641 btmp_data, n, 641 btmp_data, n,
642 pwr, &dummy_work, lwork, info 642 pwr, &dummy_work, lwork, info
643 F77_CHAR_ARG_LEN (1) 643 F77_CHAR_ARG_LEN (1)
644 F77_CHAR_ARG_LEN (1))); 644 F77_CHAR_ARG_LEN (1)));
645 645
646 if (info == 0) 646 if (info == 0)
649 Array<double> work (dim_vector (lwork, 1)); 649 Array<double> work (dim_vector (lwork, 1));
650 double *pwork = work.fortran_vec (); 650 double *pwork = work.fortran_vec ();
651 651
652 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 652 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
653 F77_CONST_CHAR_ARG2 ("U", 1), 653 F77_CONST_CHAR_ARG2 ("U", 1),
654 n, atmp_data, n, 654 n, atmp_data, n,
655 btmp_data, n, 655 btmp_data, n,
656 pwr, pwork, lwork, info 656 pwr, pwork, lwork, info
657 F77_CHAR_ARG_LEN (1) 657 F77_CHAR_ARG_LEN (1)
658 F77_CHAR_ARG_LEN (1))); 658 F77_CHAR_ARG_LEN (1)));
659 659
660 if (info < 0) 660 if (info < 0)
707 707
708 ComplexMatrix tmp = b; 708 ComplexMatrix tmp = b;
709 Complex*tmp_data = tmp.fortran_vec (); 709 Complex*tmp_data = tmp.fortran_vec ();
710 710
711 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), 711 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
712 n, tmp_data, n, 712 n, tmp_data, n,
713 info 713 info
714 F77_CHAR_ARG_LEN (1) 714 F77_CHAR_ARG_LEN (1)
715 F77_CHAR_ARG_LEN (1))); 715 F77_CHAR_ARG_LEN (1)));
716 716
717 if (a.is_hermitian () && b.is_hermitian () && info == 0) 717 if (a.is_hermitian () && b.is_hermitian () && info == 0)
743 Complex *dummy = 0; 743 Complex *dummy = 0;
744 octave_idx_type idummy = 1; 744 octave_idx_type idummy = 1;
745 745
746 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), 746 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
747 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 747 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
748 n, atmp_data, n, btmp_data, n, 748 n, atmp_data, n, btmp_data, n,
749 palpha, pbeta, dummy, idummy, 749 palpha, pbeta, dummy, idummy,
750 pv, n, &dummy_work, lwork, prwork, info 750 pv, n, &dummy_work, lwork, prwork, info
751 F77_CHAR_ARG_LEN (1) 751 F77_CHAR_ARG_LEN (1)
752 F77_CHAR_ARG_LEN (1))); 752 F77_CHAR_ARG_LEN (1)));
753 753
757 Array<Complex> work (dim_vector (lwork, 1)); 757 Array<Complex> work (dim_vector (lwork, 1));
758 Complex *pwork = work.fortran_vec (); 758 Complex *pwork = work.fortran_vec ();
759 759
760 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), 760 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
761 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 761 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
762 n, atmp_data, n, btmp_data, n, 762 n, atmp_data, n, btmp_data, n,
763 palpha, pbeta, dummy, idummy, 763 palpha, pbeta, dummy, idummy,
764 pv, n, pwork, lwork, prwork, info 764 pv, n, pwork, lwork, prwork, info
765 F77_CHAR_ARG_LEN (1) 765 F77_CHAR_ARG_LEN (1)
766 F77_CHAR_ARG_LEN (1))); 766 F77_CHAR_ARG_LEN (1)));
767 767
768 if (info < 0) 768 if (info < 0)
769 { 769 {
770 (*current_liboctave_error_handler) ("unrecoverable error in zggev"); 770 (*current_liboctave_error_handler) ("unrecoverable error in zggev");
771 return info; 771 return info;
772 } 772 }
826 Array<double> rwork (dim_vector (lrwork, 1)); 826 Array<double> rwork (dim_vector (lrwork, 1));
827 double *prwork = rwork.fortran_vec (); 827 double *prwork = rwork.fortran_vec ();
828 828
829 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 829 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
830 F77_CONST_CHAR_ARG2 ("U", 1), 830 F77_CONST_CHAR_ARG2 ("U", 1),
831 n, atmp_data, n, 831 n, atmp_data, n,
832 btmp_data, n, 832 btmp_data, n,
833 pwr, &dummy_work, lwork, 833 pwr, &dummy_work, lwork,
834 prwork, info 834 prwork, info
835 F77_CHAR_ARG_LEN (1) 835 F77_CHAR_ARG_LEN (1)
836 F77_CHAR_ARG_LEN (1))); 836 F77_CHAR_ARG_LEN (1)));
841 Array<Complex> work (dim_vector (lwork, 1)); 841 Array<Complex> work (dim_vector (lwork, 1));
842 Complex *pwork = work.fortran_vec (); 842 Complex *pwork = work.fortran_vec ();
843 843
844 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), 844 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
845 F77_CONST_CHAR_ARG2 ("U", 1), 845 F77_CONST_CHAR_ARG2 ("U", 1),
846 n, atmp_data, n, 846 n, atmp_data, n,
847 btmp_data, n, 847 btmp_data, n,
848 pwr, pwork, lwork, prwork, info 848 pwr, pwork, lwork, prwork, info
849 F77_CHAR_ARG_LEN (1) 849 F77_CHAR_ARG_LEN (1)
850 F77_CHAR_ARG_LEN (1))); 850 F77_CHAR_ARG_LEN (1)));
851 851
852 if (info < 0) 852 if (info < 0)