Mercurial > hg > octave-nkf
diff liboctave/dbleSVD.cc @ 1930:d20ab06301e8
[project @ 1996-02-11 22:30:18 by jwe]
author | jwe |
---|---|
date | Sun, 11 Feb 1996 22:30:39 +0000 |
parents | 1281a23a34dd |
children | c91f81d5f72c |
line wrap: on
line diff
--- a/liboctave/dbleSVD.cc +++ b/liboctave/dbleSVD.cc @@ -78,13 +78,14 @@ int m = a.rows (); int n = a.cols (); - double *tmp_data = dup (a.data (), a.length ()); + Matrix atmp = a; + double *tmp_data = atmp.fortran_vec (); int min_mn = m < n ? m : n; int max_mn = m > n ? m : n; - char *jobu = "A"; - char *jobv = "A"; + char jobu = 'A'; + char jobv = 'A'; int ncol_u = m; int nrow_vt = n; @@ -94,12 +95,12 @@ switch (svd_type) { case SVD::economy: - jobu = jobv ="S"; + jobu = jobv = 'S'; ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; break; case SVD::sigma_only: - jobu = jobv ="N"; + jobu = jobv = 'N'; ncol_u = nrow_vt = 1; break; @@ -109,33 +110,38 @@ type_computed = svd_type; - double *u = (*jobu == 'N' ? 0 : new double[m * ncol_u]); - double *s_vec = new double[min_mn]; - double *vt = (*jobv == 'N' ? 0 : new double[nrow_vt * n]); + if (jobu != 'N') + left_sm.resize (n, ncol_u); + + double *u = left_sm.fortran_vec (); + + sigma.resize (nrow_s, ncol_s); + double *s_vec = sigma.fortran_vec (); + + if (jobv != 'N') + right_sm.resize (nrow_vt, n); + + double *vt = right_sm.fortran_vec (); int tmp1 = 3*min_mn + max_mn; int tmp2 = 5*min_mn - 4; int lwork = tmp1 > tmp2 ? tmp1 : tmp2; - double *work = new double[lwork]; + + Array<double> work (lwork); + double *pwork = work.fortran_vec (); - F77_FCN (dgesvd, DGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u, - m, vt, nrow_vt, work, lwork, info, 1L, - 1L); - - if (*jobu != 'N') - left_sm = Matrix (u, m, ncol_u); + F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, + u, m, vt, nrow_vt, pwork, lwork, info, + 1L, 1L)); - sigma = DiagMatrix (s_vec, nrow_s, ncol_s); - - if (*jobv != 'N') + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); + else { - Matrix vt_m (vt, nrow_vt, n); - right_sm = vt_m.transpose (); + if (jobv != 'N') + right_sm = right_sm.transpose (); } - delete [] tmp_data; - delete [] work; - return info; }