Mercurial > octave-nkf

comparison liboctave/numeric/dbleCHOL.cc @ 20497:5ce959c55cc0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Propagate 'lower' in chol(a, 'lower') to underlying library function. * chol.cc (chol): Send 'L' parameter correctly when chol is called with 'lower'. * floatCHOL.cc (init): Propagate 'lower' to underlying library function. * floatCHOL.h: Modify the prototype of methods. * fMatrix.cc (inverse): Invoke chol with additional parameter. * dbleCHOL.cc (init): Propagate 'lower' to underlying library function. * dbleCHOL.h: Modify the prototype of methods. * dMatrix.cc (inverse): Invoke chol with additional parameter. * CmplxCHOL.cc (init): Propagate 'lower' to underlying library function. * CmplxCHOL.h: Modify the prototype of methods. * CMatrix.cc (inverse): Invoke chol with additional parameter.
author PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
date Sun, 24 Aug 2014 19:35:06 +0530
parents a9574e3c6e9e
children dcfbf4c1c3c8
comparison
equal deleted inserted replaced
20496:0fe7133da8ce 20497:5ce959c55cc0
85 const octave_idx_type&, double*); 85 const octave_idx_type&, double*);
86 #endif 86 #endif
87 } 87 }
88 88
89 octave_idx_type 89 octave_idx_type
90 CHOL::init (const Matrix& a, bool calc_cond) 90 CHOL::init (const Matrix& a, bool upper, bool calc_cond)
91 { 91 {
92 octave_idx_type a_nr = a.rows (); 92 octave_idx_type a_nr = a.rows ();
93 octave_idx_type a_nc = a.cols (); 93 octave_idx_type a_nc = a.cols ();
94 94
95 if (a_nr != a_nc) 95 if (a_nr != a_nc)
99 } 99 }
100 100
101 octave_idx_type n = a_nc; 101 octave_idx_type n = a_nc;
102 octave_idx_type info; 102 octave_idx_type info;
103 103
104 is_upper = upper;
105
104 chol_mat.clear (n, n); 106 chol_mat.clear (n, n);
105 for (octave_idx_type j = 0; j < n; j++) 107 if (is_upper)
106 { 108 {
107 for (octave_idx_type i = 0; i <= j; i++) 109 for (octave_idx_type j = 0; j < n; j++)
108 chol_mat.xelem (i, j) = a(i, j); 110 {
109 for (octave_idx_type i = j+1; i < n; i++) 111 for (octave_idx_type i = 0; i <= j; i++)
110 chol_mat.xelem (i, j) = 0.0; 112 chol_mat.xelem (i, j) = a (i, j);
113 for (octave_idx_type i = j + 1; i < n; i++)
114 chol_mat.xelem (i, j) = 0.0;
115 }
116 }
117 else
118 {
119 for (octave_idx_type j = 0; j < n; j++)
120 {
121 for (octave_idx_type i = 0; i < j; i++)
122 chol_mat.xelem (i, j) = 0.0;
123 for (octave_idx_type i = j; i < n; i++)
124 chol_mat.xelem (i, j) = a (i, j);
125 }
111 } 126 }
112 double *h = chol_mat.fortran_vec (); 127 double *h = chol_mat.fortran_vec ();
113 128
114 // Calculate the norm of the matrix, for later use. 129 // Calculate the norm of the matrix, for later use.
115 double anorm = 0; 130 double anorm = 0;
116 if (calc_cond) 131 if (calc_cond)
117 anorm = xnorm (a, 1); 132 anorm = xnorm (a, 1);
118 133
119 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), 134 if (is_upper)
120 n, h, n, info 135 {
121 F77_CHAR_ARG_LEN (1))); 136 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
137 n, h, n, info
138 F77_CHAR_ARG_LEN (1)));
139 }
140 else
141 {
142 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
143 n, h, n, info
144 F77_CHAR_ARG_LEN (1)));
145 }
122 146
123 xrcond = 0.0; 147 xrcond = 0.0;
124 if (info > 0) 148 if (info > 0)
125 chol_mat.resize (info - 1, info - 1); 149 chol_mat.resize (info - 1, info - 1);
126 else if (calc_cond) 150 else if (calc_cond)
130 // Now calculate the condition number for non-singular matrix. 154 // Now calculate the condition number for non-singular matrix.
131 Array<double> z (dim_vector (3*n, 1)); 155 Array<double> z (dim_vector (3*n, 1));
132 double *pz = z.fortran_vec (); 156 double *pz = z.fortran_vec ();
133 Array<octave_idx_type> iz (dim_vector (n, 1)); 157 Array<octave_idx_type> iz (dim_vector (n, 1));
134 octave_idx_type *piz = iz.fortran_vec (); 158 octave_idx_type *piz = iz.fortran_vec ();
135 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, 159 if (is_upper)
136 n, anorm, xrcond, pz, piz, dpocon_info 160 {
137 F77_CHAR_ARG_LEN (1))); 161 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
162 n, anorm, xrcond, pz, piz, dpocon_info
163 F77_CHAR_ARG_LEN (1)));
164 }
165 else
166 {
167 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
168 n, anorm, xrcond, pz, piz, dpocon_info
169 F77_CHAR_ARG_LEN (1)));
170 }
171
138 172
139 if (dpocon_info != 0) 173 if (dpocon_info != 0)
140 info = -1; 174 info = -1;
141 } 175 }
142 176
143 return info; 177 return info;
144 } 178 }
145 179
146 static Matrix 180 static Matrix
147 chol2inv_internal (const Matrix& r) 181 chol2inv_internal (const Matrix& r, bool is_upper = true)
148 { 182 {
149 Matrix retval; 183 Matrix retval;
150 184
151 octave_idx_type r_nr = r.rows (); 185 octave_idx_type r_nr = r.rows ();
152 octave_idx_type r_nc = r.cols (); 186 octave_idx_type r_nc = r.cols ();
159 Matrix tmp = r; 193 Matrix tmp = r;
160 double *v = tmp.fortran_vec (); 194 double *v = tmp.fortran_vec ();
161 195
162 if (info == 0) 196 if (info == 0)
163 { 197 {
164 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, 198 if (is_upper)
165 v, n, info 199 {
166 F77_CHAR_ARG_LEN (1))); 200 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
201 v, n, info
202 F77_CHAR_ARG_LEN (1)));
203 }
204 else
205 {
206 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
207 v, n, info
208 F77_CHAR_ARG_LEN (1)));
209 }
167 210
168 // If someone thinks of a more graceful way of doing this (or 211 // If someone thinks of a more graceful way of doing this (or
169 // faster for that matter :-)), please let me know! 212 // faster for that matter :-)), please let me know!
170 213
171 if (n > 1) 214 if (n > 1)
172 for (octave_idx_type j = 0; j < r_nc; j++) 215 {
173 for (octave_idx_type i = j+1; i < r_nr; i++) 216 if (is_upper)
174 tmp.xelem (i, j) = tmp.xelem (j, i); 217 {
218 for (octave_idx_type j = 0; j < r_nc; j++)
219 for (octave_idx_type i = j+1; i < r_nr; i++)
220 tmp.xelem (i, j) = tmp.xelem (j, i);
221 }
222 else
223 {
224 for (octave_idx_type j = 0; j < r_nc; j++)
225 for (octave_idx_type i = j+1; i < r_nr; i++)
226 tmp.xelem (j, i) = tmp.xelem (i, j);
227 }
228 }
175 229
176 retval = tmp; 230 retval = tmp;
177 } 231 }
178 } 232 }
179 else 233 else
184 238
185 // Compute the inverse of a matrix using the Cholesky factorization. 239 // Compute the inverse of a matrix using the Cholesky factorization.
186 Matrix 240 Matrix
187 CHOL::inverse (void) const 241 CHOL::inverse (void) const
188 { 242 {
189 return chol2inv_internal (chol_mat); 243 return chol2inv_internal (chol_mat, is_upper);
190 } 244 }
191 245
192 void 246 void
193 CHOL::set (const Matrix& R) 247 CHOL::set (const Matrix& R)
194 { 248 {