Mercurial > octave

comparison liboctave/dbleSVD.cc @ 1251:97eac19837dc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-04-11 15:58:32 by jwe]
author jwe
date Tue, 11 Apr 1995 15:58:32 +0000
parents b6360f2d4fa6
children bb67a902760b
comparison
equal deleted inserted replaced
1250:5cca5ae20299 1251:97eac19837dc
29 #include "mx-inlines.cc" 29 #include "mx-inlines.cc"
30 #include "f77-uscore.h" 30 #include "f77-uscore.h"
31 31
32 extern "C" 32 extern "C"
33 { 33 {
34 int F77_FCN (dgesvd) (const char*, const char*, const int*, 34 int F77_FCN (dgesvd) (const char*, const char*, const int&,
35 const int*, double*, const int*, double*, 35 const int&, double*, const int&, double*,
36 double*, const int*, double*, const int*, 36 double*, const int&, double*, const int&,
37 double*, const int*, int*, long, long); 37 double*, const int&, int&, long, long);
38 } 38 }
39 39
40 int 40 int
41 SVD::init (const Matrix& a, SVD::type svd_type) 41 SVD::init (const Matrix& a, SVD::type svd_type)
42 { 42 {
48 double *tmp_data = dup (a.data (), a.length ()); 48 double *tmp_data = dup (a.data (), a.length ());
49 49
50 int min_mn = m < n ? m : n; 50 int min_mn = m < n ? m : n;
51 int max_mn = m > n ? m : n; 51 int max_mn = m > n ? m : n;
52 52
53 char jobu = 'A'; 53 char *jobu = "A";
54 char jobv = 'A'; 54 char *jobv = "A";
55 55
56 int ncol_u = m; 56 int ncol_u = m;
57 int nrow_vt = n; 57 int nrow_vt = n;
58 int nrow_s = m; 58 int nrow_s = m;
59 int ncol_s = n; 59 int ncol_s = n;
60 60
61 if (svd_type == SVD::economy) 61 if (svd_type == SVD::economy)
62 { 62 {
63 jobu = jobv = 'S'; 63 jobu = jobv ="S";
64 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; 64 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
65 } 65 }
66 66
67 double *u = new double[m * ncol_u]; 67 double *u = new double[m * ncol_u];
68 double *s_vec = new double[min_mn]; 68 double *s_vec = new double[min_mn];
71 int tmp1 = 3*min_mn + max_mn; 71 int tmp1 = 3*min_mn + max_mn;
72 int tmp2 = 5*min_mn - 4; 72 int tmp2 = 5*min_mn - 4;
73 int lwork = tmp1 > tmp2 ? tmp1 : tmp2; 73 int lwork = tmp1 > tmp2 ? tmp1 : tmp2;
74 double *work = new double[lwork]; 74 double *work = new double[lwork];
75 75
76 F77_FCN (dgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m, 76 F77_FCN (dgesvd) (jobu, jobv, m, n, tmp_data, m, s_vec, u, m,
77 vt, &nrow_vt, work, &lwork, &info, 1L, 1L); 77 vt, nrow_vt, work, lwork, info, 1L, 1L);
78 78
79 left_sm = Matrix (u, m, ncol_u); 79 left_sm = Matrix (u, m, ncol_u);
80 sigma = DiagMatrix (s_vec, nrow_s, ncol_s); 80 sigma = DiagMatrix (s_vec, nrow_s, ncol_s);
81 Matrix vt_m (vt, nrow_vt, n); 81 Matrix vt_m (vt, nrow_vt, n);
82 right_sm = vt_m.transpose (); 82 right_sm = vt_m.transpose ();