Mercurial > hg > octave-lojdl
changeset 1930:d20ab06301e8
[project @ 1996-02-11 22:30:18 by jwe]
author | jwe |
---|---|
date | Sun, 11 Feb 1996 22:30:39 +0000 |
parents | 908f5b6676d7 |
children | c91f81d5f72c |
files | liboctave/CmplxSCHUR.cc liboctave/CmplxSCHUR.h liboctave/CmplxSVD.cc liboctave/CmplxSVD.h liboctave/dbleSCHUR.cc liboctave/dbleSCHUR.h liboctave/dbleSVD.cc liboctave/dbleSVD.h |
diffstat | 8 files changed, 81 insertions(+), 58 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxSCHUR.cc +++ b/liboctave/CmplxSCHUR.cc @@ -70,20 +70,21 @@ return -1; } - char *jobvs = "V"; - char *sense = "N"; - char *sort = "N"; + char jobvs = 'V'; + char sense = 'N'; + char sort = 'N'; char ord_char = ord.empty () ? 'U' : ord[0]; if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = "S"; + sort = 'S'; - select_function selector= 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; + else + selector = 0; int n = a_nc; int lwork = 8 * n; @@ -116,7 +117,7 @@ int *pbwork = bwork.fortran_vec (); - F77_XFCN (zgeesx, ZGEESX, (jobvs, sort, selector, sense, n, s, n, + F77_XFCN (zgeesx, ZGEESX, (&jobvs, &sort, selector, &sense, n, s, n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork, info, 1L, 1L));
--- a/liboctave/CmplxSCHUR.h +++ b/liboctave/CmplxSCHUR.h @@ -82,6 +82,8 @@ ComplexMatrix schur_mat; ComplexMatrix unitary_mat; + select_function selector; + int init (const ComplexMatrix& a, const string& ord); };
--- a/liboctave/CmplxSVD.cc +++ b/liboctave/CmplxSVD.cc @@ -83,8 +83,8 @@ 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 +94,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,34 +109,41 @@ type_computed = svd_type; - Complex *u = (*jobu == 'N' ? 0 : new Complex[m * ncol_u]); - double *s_vec = new double[min_mn]; - Complex *vt = (*jobv == 'N' ? 0 : new Complex[nrow_vt * n]); + if (jobu != 'N') + left_sm.resize (m, ncol_u); + + Complex *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); + + Complex *vt = right_sm.fortran_vec (); int lwork = 2*min_mn + max_mn; - Complex *work = new Complex[lwork]; + + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); int lrwork = 5*max_mn; - double *rwork = new double[lrwork]; + + Array<double> rwork (lrwork); + double *prwork = rwork.fortran_vec (); - F77_FCN (zgesvd, ZGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u, - m, vt, nrow_vt, work, lwork, rwork, info, - 1L, 1L); - - if (*jobu != 'N') - left_sm = ComplexMatrix (u, m, ncol_u); + F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, u, + m, vt, nrow_vt, pwork, lwork, prwork, 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 zgesvd"); + else { - ComplexMatrix vt_m (vt, nrow_vt, n); - right_sm = vt_m.hermitian (); + if (jobv != 'N') + right_sm = right_sm.hermitian (); } - delete [] tmp_data; - delete [] work; - return info; }
--- a/liboctave/CmplxSVD.h +++ b/liboctave/CmplxSVD.h @@ -66,6 +66,8 @@ return *this; } + ~ComplexSVD (void) { } + DiagMatrix singular_values (void) const { return sigma; } ComplexMatrix left_singular_matrix (void) const;
--- a/liboctave/dbleSCHUR.cc +++ b/liboctave/dbleSCHUR.cc @@ -71,20 +71,21 @@ return -1; } - char *jobvs = "V"; - char *sense = "N"; - char *sort = "N"; + char jobvs = 'V'; + char sense = 'N'; + char sort = 'N'; char ord_char = ord.empty () ? 'U' : ord[0]; if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = "S"; + sort = 'S'; - select_function selector = 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; + else + selector = 0; int n = a_nc; int lwork = 8 * n; @@ -124,7 +125,7 @@ int *piwork = iwork.fortran_vec (); - F77_XFCN (dgeesx, DGEESX, (jobvs, sort, selector, sense, n, s, + F77_XFCN (dgeesx, DGEESX, (&jobvs, &sort, selector, &sense, n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork, liwork, pbwork, info, 1L, 1L));
--- a/liboctave/dbleSCHUR.h +++ b/liboctave/dbleSCHUR.h @@ -82,6 +82,8 @@ Matrix schur_mat; Matrix unitary_mat; + select_function selector; + int init (const Matrix& a, const string& ord); };
--- 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; }