Mercurial > octave-nkf

comparison liboctave/CMatrix.cc @ 4552:6f3382e08a52

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2003-10-27 20:38:02 by jwe]
author jwe
date Mon, 27 Oct 2003 20:38:03 +0000
parents 508238e65af7
children 7b957b442818
comparison
equal deleted inserted replaced
4551:2c619e5138fd 4552:6f3382e08a52
61 61
62 // Fortran functions we call. 62 // Fortran functions we call.
63 63
64 extern "C" 64 extern "C"
65 { 65 {
66 int F77_FUNC (zgebal, ZGEBAL) (const char*, const int&, Complex*, 66 F77_RET_T
67 const int&, int&, int&, double*, int&, 67 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
68 long, long); 68 const int&, Complex*, const int&, int&,
69 69 int&, double*, int&
70 int F77_FUNC (dgebak, DGEBAK) (const char*, const char*, const int&, 70 F77_CHAR_ARG_LEN_DECL);
71 const int&, const int&, double*, 71
72 const int&, double*, const int&, 72 F77_RET_T
73 int&, long, long); 73 F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
74 74 F77_CONST_CHAR_ARG_DECL,
75 int F77_FUNC (zgemm, ZGEMM) (const char*, const char*, const int&, 75 const int&, const int&, const int&, double*,
76 const int&, const int&, const Complex&, 76 const int&, double*, const int&, int&
77 const Complex*, const int&, 77 F77_CHAR_ARG_LEN_DECL
78 const Complex*, const int&, 78 F77_CHAR_ARG_LEN_DECL);
79 const Complex&, Complex*, const int&, 79
80 long, long); 80 F77_RET_T
81 81 F77_FUNC (zgemm, ZGEMM) (F77_CONST_CHAR_ARG_DECL,
82 int F77_FUNC (zgetrf, ZGETRF) (const int&, const int&, Complex*, const int&, 82 F77_CONST_CHAR_ARG_DECL,
83 int*, int&); 83 const int&, const int&, const int&,
84 84 const Complex&, const Complex*, const int&,
85 int F77_FUNC (zgetrs, ZGETRS) (const char*, const int&, const int&, 85 const Complex*, const int&, const Complex&,
86 Complex*, const int&, 86 Complex*, const int&
87 const int*, Complex*, const int&, int&); 87 F77_CHAR_ARG_LEN_DECL
88 88 F77_CHAR_ARG_LEN_DECL);
89 int F77_FUNC (zgetri, ZGETRI) (const int&, Complex*, const int&, const int*, 89
90 Complex*, const int&, int&); 90 F77_RET_T
91 91 F77_FUNC (zgetrf, ZGETRF) (const int&, const int&, Complex*, const int&,
92 int F77_FUNC (zgecon, ZGECON) (const char*, const int&, Complex*, 92 int*, int&);
93 const int&, const double&, double&, 93
94 Complex*, double*, int&); 94 F77_RET_T
95 95 F77_FUNC (zgetrs, ZGETRS) (F77_CONST_CHAR_ARG_DECL,
96 int F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&, 96 const int&, const int&, Complex*, const int&,
97 Complex*, const int&, Complex*, 97 const int*, Complex*, const int&, int&
98 const int&, double*, double&, int&, 98 F77_CHAR_ARG_LEN_DECL);
99 Complex*, const int&, double*, int&); 99
100 F77_RET_T
101 F77_FUNC (zgetri, ZGETRI) (const int&, Complex*, const int&, const int*,
102 Complex*, const int&, int&);
103
104 F77_RET_T
105 F77_FUNC (zgecon, ZGECON) (F77_CONST_CHAR_ARG_DECL,
106 const int&, Complex*,
107 const int&, const double&, double&,
108 Complex*, double*, int&
109 F77_CHAR_ARG_LEN_DECL);
110
111 F77_RET_T
112 F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&,
113 Complex*, const int&, Complex*,
114 const int&, double*, double&, int&,
115 Complex*, const int&, double*, int&);
100 116
101 // Note that the original complex fft routines were not written for 117 // Note that the original complex fft routines were not written for
102 // double complex arguments. They have been modified by adding an 118 // double complex arguments. They have been modified by adding an
103 // implicit double precision (a-h,o-z) statement at the beginning of 119 // implicit double precision (a-h,o-z) statement at the beginning of
104 // each subroutine. 120 // each subroutine.
105 121
106 int F77_FUNC (cffti, CFFTI) (const int&, Complex*); 122 F77_RET_T
107 123 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
108 int F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*); 124
109 125 F77_RET_T
110 int F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*); 126 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
111 127
112 int F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&, 128 F77_RET_T
113 double&, Complex&, Complex&); 129 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
114 130
115 int F77_FUNC (ztrsyl, ZTRSYL) (const char*, const char*, const int&, 131 F77_RET_T
116 const int&, const int&, 132 F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&,
117 const Complex*, const int&, 133 double&, Complex&, Complex&);
118 const Complex*, const int&, 134
119 const Complex*, const int&, double&, 135 F77_RET_T
120 int&, long, long); 136 F77_FUNC (ztrsyl, ZTRSYL) (F77_CONST_CHAR_ARG_DECL,
121 137 F77_CONST_CHAR_ARG_DECL,
122 int F77_FUNC (xzlange, XZLANGE) (const char*, const int&, 138 const int&, const int&, const int&,
123 const int&, const Complex*, 139 const Complex*, const int&,
124 const int&, double*, double&); 140 const Complex*, const int&,
141 const Complex*, const int&, double&, int&
142 F77_CHAR_ARG_LEN_DECL
143 F77_CHAR_ARG_LEN_DECL);
144
145 F77_RET_T
146 F77_FUNC (xzlange, XZLANGE) (F77_CONST_CHAR_ARG_DECL,
147 const int&, const int&, const Complex*,
148 const int&, double*, double&
149 F77_CHAR_ARG_LEN_DECL);
125 } 150 }
126 151
127 static const Complex Complex_NaN_result (octave_NaN, octave_NaN); 152 static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
128 153
129 // Complex Matrix class 154 // Complex Matrix class
1000 { 1025 {
1001 // Now calculate the condition number for non-singular matrix. 1026 // Now calculate the condition number for non-singular matrix.
1002 char job = '1'; 1027 char job = '1';
1003 Array<double> rz (2 * nc); 1028 Array<double> rz (2 * nc);
1004 double *prz = rz.fortran_vec (); 1029 double *prz = rz.fortran_vec ();
1005 F77_XFCN (zgecon, ZGECON, (&job, nc, tmp_data, nr, anorm, 1030 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1006 rcond, pz, prz, info)); 1031 nc, tmp_data, nr, anorm,
1032 rcond, pz, prz, info
1033 F77_CHAR_ARG_LEN (1)));
1007 1034
1008 if (f77_exception_encountered) 1035 if (f77_exception_encountered)
1009 (*current_liboctave_error_handler) 1036 (*current_liboctave_error_handler)
1010 ("unrecoverable error in zgecon"); 1037 ("unrecoverable error in zgecon");
1011 1038
1016 if (info == -1 && ! force) 1043 if (info == -1 && ! force)
1017 retval = *this; // Restore contents. 1044 retval = *this; // Restore contents.
1018 else 1045 else
1019 { 1046 {
1020 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, 1047 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
1021 pz, lwork, info)); 1048 pz, lwork, info));
1022 1049
1023 if (f77_exception_encountered) 1050 if (f77_exception_encountered)
1024 (*current_liboctave_error_handler) 1051 (*current_liboctave_error_handler)
1025 ("unrecoverable error in zgetri"); 1052 ("unrecoverable error in zgetri");
1026 1053
1463 Array<Complex> z (2*nr); 1490 Array<Complex> z (2*nr);
1464 Complex *pz = z.fortran_vec (); 1491 Complex *pz = z.fortran_vec ();
1465 Array<double> rz (2*nr); 1492 Array<double> rz (2*nr);
1466 double *prz = rz.fortran_vec (); 1493 double *prz = rz.fortran_vec ();
1467 1494
1468 F77_XFCN (zgecon, ZGECON, (&job, nc, tmp_data, nr, anorm, 1495 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1469 rcond, pz, prz, info)); 1496 nc, tmp_data, nr, anorm,
1497 rcond, pz, prz, info
1498 F77_CHAR_ARG_LEN (1)));
1470 1499
1471 if (f77_exception_encountered) 1500 if (f77_exception_encountered)
1472 (*current_liboctave_error_handler) 1501 (*current_liboctave_error_handler)
1473 ("unrecoverable error in zgecon"); 1502 ("unrecoverable error in zgecon");
1474 } 1503 }
1607 } 1636 }
1608 else 1637 else
1609 { 1638 {
1610 // Now calculate the condition number for non-singular matrix. 1639 // Now calculate the condition number for non-singular matrix.
1611 char job = '1'; 1640 char job = '1';
1612 F77_XFCN (zgecon, ZGECON, (&job, nc, tmp_data, nr, anorm, 1641 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1613 rcond, pz, prz, info)); 1642 nc, tmp_data, nr, anorm,
1643 rcond, pz, prz, info
1644 F77_CHAR_ARG_LEN (1)));
1614 1645
1615 if (f77_exception_encountered) 1646 if (f77_exception_encountered)
1616 (*current_liboctave_error_handler) 1647 (*current_liboctave_error_handler)
1617 ("unrecoverable error in zgecon"); 1648 ("unrecoverable error in zgecon");
1618 1649
1638 Complex *result = retval.fortran_vec (); 1669 Complex *result = retval.fortran_vec ();
1639 1670
1640 int b_nc = b.cols (); 1671 int b_nc = b.cols ();
1641 1672
1642 char job = 'N'; 1673 char job = 'N';
1643 F77_XFCN (zgetrs, ZGETRS, (&job, nr, b_nc, tmp_data, nr, 1674 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1644 pipvt, result, b.rows(), info)); 1675 nr, b_nc, tmp_data, nr,
1676 pipvt, result, b.rows(), info
1677 F77_CHAR_ARG_LEN (1)));
1645 1678
1646 if (f77_exception_encountered) 1679 if (f77_exception_encountered)
1647 (*current_liboctave_error_handler) 1680 (*current_liboctave_error_handler)
1648 ("unrecoverable error in zgetrs"); 1681 ("unrecoverable error in zgetrs");
1649 } 1682 }
1756 } 1789 }
1757 else 1790 else
1758 { 1791 {
1759 // Now calculate the condition number for non-singular matrix. 1792 // Now calculate the condition number for non-singular matrix.
1760 char job = '1'; 1793 char job = '1';
1761 F77_XFCN (zgecon, ZGECON, (&job, nc, tmp_data, nr, anorm, 1794 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1762 rcond, pz, prz, info)); 1795 nc, tmp_data, nr, anorm,
1796 rcond, pz, prz, info
1797 F77_CHAR_ARG_LEN (1)));
1763 1798
1764 if (f77_exception_encountered) 1799 if (f77_exception_encountered)
1765 (*current_liboctave_error_handler) 1800 (*current_liboctave_error_handler)
1766 ("unrecoverable error in zgecon"); 1801 ("unrecoverable error in zgecon");
1767 1802
1785 { 1820 {
1786 retval = b; 1821 retval = b;
1787 Complex *result = retval.fortran_vec (); 1822 Complex *result = retval.fortran_vec ();
1788 1823
1789 char job = 'N'; 1824 char job = 'N';
1790 F77_XFCN (zgetrs, ZGETRS, (&job, nr, 1, tmp_data, nr, pipvt, 1825 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1791 result, b.length(), info)); 1826 nr, 1, tmp_data, nr, pipvt,
1827 result, b.length(), info
1828 F77_CHAR_ARG_LEN (1)));
1792 1829
1793 if (f77_exception_encountered) 1830 if (f77_exception_encountered)
1794 (*current_liboctave_error_handler) 1831 (*current_liboctave_error_handler)
1795 ("unrecoverable error in zgetrs"); 1832 ("unrecoverable error in zgetrs");
1796 1833
2077 2114
2078 // XXX FIXME XXX -- should pass job as a parameter in expm 2115 // XXX FIXME XXX -- should pass job as a parameter in expm
2079 2116
2080 // Permute first 2117 // Permute first
2081 char job = 'P'; 2118 char job = 'P';
2082 F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi, 2119 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2083 dpermute.fortran_vec (), info, 1L, 1L)); 2120 nc, mp, nc, ilo, ihi,
2121 dpermute.fortran_vec (), info
2122 F77_CHAR_ARG_LEN (1)));
2084 2123
2085 if (f77_exception_encountered) 2124 if (f77_exception_encountered)
2086 { 2125 {
2087 (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); 2126 (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
2088 return retval; 2127 return retval;
2089 } 2128 }
2090 2129
2091 // then scale 2130 // then scale
2092 job = 'S'; 2131 job = 'S';
2093 F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilos, ihis, 2132 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2094 dscale.fortran_vec (), info, 1L, 1L)); 2133 nc, mp, nc, ilos, ihis,
2134 dscale.fortran_vec (), info
2135 F77_CHAR_ARG_LEN (1)));
2095 2136
2096 if (f77_exception_encountered) 2137 if (f77_exception_encountered)
2097 { 2138 {
2098 (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); 2139 (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
2099 return retval; 2140 return retval;
2102 // Preconditioning step 3: scaling. 2143 // Preconditioning step 3: scaling.
2103 2144
2104 ColumnVector work (nc); 2145 ColumnVector work (nc);
2105 double inf_norm; 2146 double inf_norm;
2106 2147
2107 F77_XFCN (xzlange, XZLANGE, ("I", nc, nc, m.fortran_vec (), nc, 2148 F77_XFCN (xzlange, XZLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
2108 work.fortran_vec (), inf_norm)); 2149 nc, nc, m.fortran_vec (), nc,
2150 work.fortran_vec (), inf_norm
2151 F77_CHAR_ARG_LEN (1)));
2109 2152
2110 if (f77_exception_encountered) 2153 if (f77_exception_encountered)
2111 { 2154 {
2112 (*current_liboctave_error_handler) ("unrecoverable error in zlange"); 2155 (*current_liboctave_error_handler) ("unrecoverable error in zlange");
2113 return retval; 2156 return retval;
2244 int a_len = a.length (); 2287 int a_len = a.length ();
2245 2288
2246 retval.resize (len, a_len); 2289 retval.resize (len, a_len);
2247 Complex *c = retval.fortran_vec (); 2290 Complex *c = retval.fortran_vec ();
2248 2291
2249 F77_XFCN (zgemm, ZGEMM, ("N", "N", len, a_len, 1, 1.0, 2292 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
2250 v.data (), len, a.data (), 1, 0.0, 2293 F77_CONST_CHAR_ARG2 ("N", 1),
2251 c, len, 1L, 1L)); 2294 len, a_len, 1, 1.0, v.data (), len,
2295 a.data (), 1, 0.0, c, len
2296 F77_CHAR_ARG_LEN (1)
2297 F77_CHAR_ARG_LEN (1)));
2252 2298
2253 if (f77_exception_encountered) 2299 if (f77_exception_encountered)
2254 (*current_liboctave_error_handler) 2300 (*current_liboctave_error_handler)
2255 ("unrecoverable error in zgemm"); 2301 ("unrecoverable error in zgemm");
2256 } 2302 }
3128 3174
3129 Complex *pa = sch_a.fortran_vec (); 3175 Complex *pa = sch_a.fortran_vec ();
3130 Complex *pb = sch_b.fortran_vec (); 3176 Complex *pb = sch_b.fortran_vec ();
3131 Complex *px = cx.fortran_vec (); 3177 Complex *px = cx.fortran_vec ();
3132 3178
3133 F77_XFCN (ztrsyl, ZTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, 3179 F77_XFCN (ztrsyl, ZTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
3134 b_nr, px, a_nr, scale, 3180 F77_CONST_CHAR_ARG2 ("N", 1),
3135 info, 1L, 1L)); 3181 1, a_nr, b_nr, pa, a_nr, pb,
3182 b_nr, px, a_nr, scale, info
3183 F77_CHAR_ARG_LEN (1)
3184 F77_CHAR_ARG_LEN (1)));
3136 3185
3137 if (f77_exception_encountered) 3186 if (f77_exception_encountered)
3138 (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl"); 3187 (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl");
3139 else 3188 else
3140 { 3189 {
3183 int lda = a.rows (); 3232 int lda = a.rows ();
3184 3233
3185 retval.resize (nr, a_nc); 3234 retval.resize (nr, a_nc);
3186 Complex *c = retval.fortran_vec (); 3235 Complex *c = retval.fortran_vec ();
3187 3236
3188 F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0, 3237 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
3189 m.data (), ld, a.data (), lda, 0.0, 3238 F77_CONST_CHAR_ARG2 ("N", 1),
3190 c, nr, 1L, 1L)); 3239 nr, a_nc, nc, 1.0, m.data (),
3240 ld, a.data (), lda, 0.0, c, nr
3241 F77_CHAR_ARG_LEN (1)
3242 F77_CHAR_ARG_LEN (1)));
3191 3243
3192 if (f77_exception_encountered) 3244 if (f77_exception_encountered)
3193 (*current_liboctave_error_handler) 3245 (*current_liboctave_error_handler)
3194 ("unrecoverable error in zgemm"); 3246 ("unrecoverable error in zgemm");
3195 } 3247 }