Mercurial > hg > octave-nkf

comparison liboctave/fCmplxSVD.cc @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children eb63fbe60fab
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 /*
2
3 Copyright (C) 1994, 1995, 1996, 1997, 1999, 2002, 2003, 2004, 2005,
4 2007 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27
28 #include "fCmplxSVD.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31
32 extern "C"
33 {
34 F77_RET_T
35 F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
36 F77_CONST_CHAR_ARG_DECL,
37 const octave_idx_type&, const octave_idx_type&, FloatComplex*,
38 const octave_idx_type&, float*, FloatComplex*, const octave_idx_type&,
39 FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&,
40 float*, octave_idx_type&
41 F77_CHAR_ARG_LEN_DECL
42 F77_CHAR_ARG_LEN_DECL);
43 }
44
45 FloatComplexMatrix
46 FloatComplexSVD::left_singular_matrix (void) const
47 {
48 if (type_computed == SVD::sigma_only)
49 {
50 (*current_liboctave_error_handler)
51 ("FloatComplexSVD: U not computed because type == SVD::sigma_only");
52 return FloatComplexMatrix ();
53 }
54 else
55 return left_sm;
56 }
57
58 FloatComplexMatrix
59 FloatComplexSVD::right_singular_matrix (void) const
60 {
61 if (type_computed == SVD::sigma_only)
62 {
63 (*current_liboctave_error_handler)
64 ("FloatComplexSVD: V not computed because type == SVD::sigma_only");
65 return FloatComplexMatrix ();
66 }
67 else
68 return right_sm;
69 }
70
71 octave_idx_type
72 FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type)
73 {
74 octave_idx_type info;
75
76 octave_idx_type m = a.rows ();
77 octave_idx_type n = a.cols ();
78
79 FloatComplexMatrix atmp = a;
80 FloatComplex *tmp_data = atmp.fortran_vec ();
81
82 octave_idx_type min_mn = m < n ? m : n;
83 octave_idx_type max_mn = m > n ? m : n;
84
85 char jobu = 'A';
86 char jobv = 'A';
87
88 octave_idx_type ncol_u = m;
89 octave_idx_type nrow_vt = n;
90 octave_idx_type nrow_s = m;
91 octave_idx_type ncol_s = n;
92
93 switch (svd_type)
94 {
95 case SVD::economy:
96 jobu = jobv = 'S';
97 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
98 break;
99
100 case SVD::sigma_only:
101
102 // Note: for this case, both jobu and jobv should be 'N', but
103 // there seems to be a bug in dgesvd from Lapack V2.0. To
104 // demonstrate the bug, set both jobu and jobv to 'N' and find
105 // the singular values of [eye(3), eye(3)]. The result is
106 // [-sqrt(2), -sqrt(2), -sqrt(2)].
107 //
108 // For Lapack 3.0, this problem seems to be fixed.
109
110 jobu = 'N';
111 jobv = 'N';
112 ncol_u = nrow_vt = 1;
113 break;
114
115 default:
116 break;
117 }
118
119 type_computed = svd_type;
120
121 if (! (jobu == 'N' || jobu == 'O'))
122 left_sm.resize (m, ncol_u);
123
124 FloatComplex *u = left_sm.fortran_vec ();
125
126 sigma.resize (nrow_s, ncol_s);
127 float *s_vec = sigma.fortran_vec ();
128
129 if (! (jobv == 'N' || jobv == 'O'))
130 right_sm.resize (nrow_vt, n);
131
132 FloatComplex *vt = right_sm.fortran_vec ();
133
134 octave_idx_type lrwork = 5*max_mn;
135
136 Array<float> rwork (lrwork);
137
138 // Ask ZGESVD what the dimension of WORK should be.
139
140 octave_idx_type lwork = -1;
141
142 Array<FloatComplex> work (1);
143
144 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
145 F77_CONST_CHAR_ARG2 (&jobv, 1),
146 m, n, tmp_data, m, s_vec, u, m, vt,
147 nrow_vt, work.fortran_vec (), lwork,
148 rwork.fortran_vec (), info
149 F77_CHAR_ARG_LEN (1)
150 F77_CHAR_ARG_LEN (1)));
151
152 lwork = static_cast<octave_idx_type> (work(0).real ());
153 work.resize (lwork);
154
155 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
156 F77_CONST_CHAR_ARG2 (&jobv, 1),
157 m, n, tmp_data, m, s_vec, u, m, vt,
158 nrow_vt, work.fortran_vec (), lwork,
159 rwork.fortran_vec (), info
160 F77_CHAR_ARG_LEN (1)
161 F77_CHAR_ARG_LEN (1)));
162
163 if (! (jobv == 'N' || jobv == 'O'))
164 right_sm = right_sm.hermitian ();
165
166 return info;
167 }
168
169 /*
170 ;;; Local Variables: ***
171 ;;; mode: C++ ***
172 ;;; End: ***
173 */