Mercurial > hg > octave-lyh

comparison liboctave/dbleSVD.cc @ 1545:8bdfa6fe5d69

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-10-06 06:13:28 by jwe]
author jwe
date Fri, 06 Oct 1995 06:14:45 +0000
parents f1931fc63ce9
children b0b22b6ce22f
comparison
equal deleted inserted replaced
1544:f1931fc63ce9 1545:8bdfa6fe5d69
41 const int&, double*, const int&, int&, 41 const int&, double*, const int&, int&,
42 long, long); 42 long, long);
43 } 43 }
44 44
45 Matrix 45 Matrix
46 left_singular_matrix (void) const 46 SVD::left_singular_matrix (void) const
47 { 47 {
48 if (type_computed == SVD::sigma_only) 48 if (type_computed == SVD::sigma_only)
49 { 49 {
50 (*current_liboctave_error_handler) 50 (*current_liboctave_error_handler)
51 ("ComplexSVD: U not computed because type == SVD::sigma_only"); 51 ("ComplexSVD: U not computed because type == SVD::sigma_only");
54 else 54 else
55 return left_sm; 55 return left_sm;
56 } 56 }
57 57
58 Matrix 58 Matrix
59 right_singular_matrix (void) const 59 SVD::right_singular_matrix (void) const
60 { 60 {
61 if (type_computed == SVD::sigma_only) 61 if (type_computed == SVD::sigma_only)
62 { 62 {
63 (*current_liboctave_error_handler) 63 (*current_liboctave_error_handler)
64 ("ComplexSVD: V not computed because type == SVD::sigma_only"); 64 ("ComplexSVD: V not computed because type == SVD::sigma_only");
96 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; 96 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
97 break; 97 break;
98 98
99 case SVD::sigma_only: 99 case SVD::sigma_only:
100 jobu = jobv ="N"; 100 jobu = jobv ="N";
101 ncol_u = nrow_vt = 0; 101 ncol_u = nrow_vt = 1;
102 break; 102 break;
103 103
104 default: 104 default:
105 break; 105 break;
106 } 106 }
107 107
108 type_computed = svd_type; 108 type_computed = svd_type;
109 109
110 double *u = (ncol_u > 0 ? new double[m * ncol_u] : 0); 110 double *u = (*jobu == 'N' ? 0 : new double[m * ncol_u]);
111 double *s_vec = new double[min_mn]; 111 double *s_vec = new double[min_mn];
112 double *vt = (ncol_vt > 0 ? new double[nrow_vt * n] : 0); 112 double *vt = (*jobv == 'N' ? 0 : new double[nrow_vt * n]);
113 113
114 int tmp1 = 3*min_mn + max_mn; 114 int tmp1 = 3*min_mn + max_mn;
115 int tmp2 = 5*min_mn - 4; 115 int tmp2 = 5*min_mn - 4;
116 int lwork = tmp1 > tmp2 ? tmp1 : tmp2; 116 int lwork = tmp1 > tmp2 ? tmp1 : tmp2;
117 double *work = new double[lwork]; 117 double *work = new double[lwork];
118 118
119 F77_FCN (dgesvd, DGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u, 119 F77_FCN (dgesvd, DGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u,
120 m, vt, nrow_vt, work, lwork, info, 1L, 120 m, vt, nrow_vt, work, lwork, info, 1L,
121 1L); 121 1L);
122 122
123 if (ncol_u > 0) 123 if (*jobu != 'N')
124 left_sm = Matrix (u, m, ncol_u); 124 left_sm = Matrix (u, m, ncol_u);
125 125
126 sigma = DiagMatrix (s_vec, nrow_s, ncol_s); 126 sigma = DiagMatrix (s_vec, nrow_s, ncol_s);
127 127
128 if (nrow_vt > 0) 128 if (*jobv != 'N')
129 { 129 {
130 Matrix vt_m (vt, nrow_vt, n); 130 Matrix vt_m (vt, nrow_vt, n);
131 right_sm = vt_m.transpose (); 131 right_sm = vt_m.transpose ();
132 } 132 }
133 133