Mercurial > octave-nkf

comparison liboctave/SparseQR.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents feb90c7cfc16
children 72c96de7a403
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
72 CXSPARSE_DNAME (_sfree) (S); 72 CXSPARSE_DNAME (_sfree) (S);
73 CXSPARSE_DNAME (_nfree) (N); 73 CXSPARSE_DNAME (_nfree) (N);
74 #endif 74 #endif
75 } 75 }
76 76
77 SparseMatrix 77 SparseMatrix
78 SparseQR::SparseQR_rep::V (void) const 78 SparseQR::SparseQR_rep::V (void) const
79 { 79 {
80 #ifdef HAVE_CXSPARSE 80 #ifdef HAVE_CXSPARSE
81 // Drop zeros from V and sort 81 // Drop zeros from V and sort
82 // FIXME Is the double transpose to sort necessary? 82 // FIXME Is the double transpose to sort necessary?
102 #else 102 #else
103 return SparseMatrix (); 103 return SparseMatrix ();
104 #endif 104 #endif
105 } 105 }
106 106
107 ColumnVector 107 ColumnVector
108 SparseQR::SparseQR_rep::Pinv (void) const 108 SparseQR::SparseQR_rep::Pinv (void) const
109 { 109 {
110 #ifdef HAVE_CXSPARSE 110 #ifdef HAVE_CXSPARSE
111 ColumnVector ret(N->L->m); 111 ColumnVector ret(N->L->m);
112 for (octave_idx_type i = 0; i < N->L->m; i++) 112 for (octave_idx_type i = 0; i < N->L->m; i++)
119 #else 119 #else
120 return ColumnVector (); 120 return ColumnVector ();
121 #endif 121 #endif
122 } 122 }
123 123
124 ColumnVector 124 ColumnVector
125 SparseQR::SparseQR_rep::P (void) const 125 SparseQR::SparseQR_rep::P (void) const
126 { 126 {
127 #ifdef HAVE_CXSPARSE 127 #ifdef HAVE_CXSPARSE
128 ColumnVector ret(N->L->m); 128 ColumnVector ret(N->L->m);
129 for (octave_idx_type i = 0; i < N->L->m; i++) 129 for (octave_idx_type i = 0; i < N->L->m; i++)
136 #else 136 #else
137 return ColumnVector (); 137 return ColumnVector ();
138 #endif 138 #endif
139 } 139 }
140 140
141 SparseMatrix 141 SparseMatrix
142 SparseQR::SparseQR_rep::R (const bool econ) const 142 SparseQR::SparseQR_rep::R (const bool econ) const
143 { 143 {
144 #ifdef HAVE_CXSPARSE 144 #ifdef HAVE_CXSPARSE
145 // Drop zeros from R and sort 145 // Drop zeros from R and sort
146 // FIXME Is the double transpose to sort necessary? 146 // FIXME Is the double transpose to sort necessary?
288 else if (nr == 0 || nc == 0 || b_nc == 0) 288 else if (nr == 0 || nc == 0 || b_nc == 0)
289 x = Matrix (nc, b_nc, 0.0); 289 x = Matrix (nc, b_nc, 0.0);
290 else if (nr >= nc) 290 else if (nr >= nc)
291 { 291 {
292 SparseQR q (a, 3); 292 SparseQR q (a, 3);
293 if (! q.ok ()) 293 if (! q.ok ())
294 return Matrix(); 294 return Matrix();
295 x.resize(nc, b_nc); 295 x.resize(nc, b_nc);
296 double *vec = x.fortran_vec(); 296 double *vec = x.fortran_vec();
297 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); 297 OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2);
298 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 298 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
299 i++, idx+=nc, bidx+=b_nr) 299 i++, idx+=nc, bidx+=b_nr)
300 { 300 {
301 octave_quit (); 301 octave_quit ();
302 for (octave_idx_type j = nr; j < q.S()->m2; j++) 302 for (octave_idx_type j = nr; j < q.S()->m2; j++)
303 buf[j] = 0.; 303 buf[j] = 0.;
334 return Matrix(); 334 return Matrix();
335 x.resize(nc, b_nc); 335 x.resize(nc, b_nc);
336 double *vec = x.fortran_vec(); 336 double *vec = x.fortran_vec();
337 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2); 337 volatile octave_idx_type nbuf = (nc > q.S()->m2 ? nc : q.S()->m2);
338 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 338 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
339 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 339 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
340 i++, idx+=nc, bidx+=b_nr) 340 i++, idx+=nc, bidx+=b_nr)
341 { 341 {
342 octave_quit (); 342 octave_quit ();
343 for (octave_idx_type j = nr; j < nbuf; j++) 343 for (octave_idx_type j = nr; j < nbuf; j++)
344 buf[j] = 0.; 344 buf[j] = 0.;
392 else if (nr == 0 || nc == 0 || b_nc == 0) 392 else if (nr == 0 || nc == 0 || b_nc == 0)
393 x = SparseMatrix (nc, b_nc); 393 x = SparseMatrix (nc, b_nc);
394 else if (nr >= nc) 394 else if (nr >= nc)
395 { 395 {
396 SparseQR q (a, 3); 396 SparseQR q (a, 3);
397 if (! q.ok ()) 397 if (! q.ok ())
398 return SparseMatrix(); 398 return SparseMatrix();
399 x = SparseMatrix (nc, b_nc, b.nnz()); 399 x = SparseMatrix (nc, b_nc, b.nnz());
400 x.xcidx(0) = 0; 400 x.xcidx(0) = 0;
401 x_nz = b.nnz(); 401 x_nz = b.nnz();
402 ii = 0; 402 ii = 0;
710 else if (nr == 0 || nc == 0 || b_nc == 0) 710 else if (nr == 0 || nc == 0 || b_nc == 0)
711 x = SparseComplexMatrix (nc, b_nc); 711 x = SparseComplexMatrix (nc, b_nc);
712 else if (nr >= nc) 712 else if (nr >= nc)
713 { 713 {
714 SparseQR q (a, 3); 714 SparseQR q (a, 3);
715 if (! q.ok ()) 715 if (! q.ok ())
716 return SparseComplexMatrix(); 716 return SparseComplexMatrix();
717 x = SparseComplexMatrix (nc, b_nc, b.nnz()); 717 x = SparseComplexMatrix (nc, b_nc, b.nnz());
718 x.xcidx(0) = 0; 718 x.xcidx(0) = 0;
719 x_nz = b.nnz(); 719 x_nz = b.nnz();
720 ii = 0; 720 ii = 0;
899 #else 899 #else
900 return SparseComplexMatrix (); 900 return SparseComplexMatrix ();
901 #endif 901 #endif
902 } 902 }
903 903
904 Matrix 904 Matrix
905 qrsolve(const SparseMatrix &a, const MArray<double> &b, 905 qrsolve(const SparseMatrix &a, const MArray<double> &b,
906 octave_idx_type &info) 906 octave_idx_type &info)
907 { 907 {
908 return qrsolve (a, Matrix (b), info); 908 return qrsolve (a, Matrix (b), info);
909 } 909 }
910 910
911 ComplexMatrix 911 ComplexMatrix
912 qrsolve(const SparseMatrix &a, const MArray<Complex> &b, 912 qrsolve(const SparseMatrix &a, const MArray<Complex> &b,
913 octave_idx_type &info) 913 octave_idx_type &info)
914 { 914 {
915 return qrsolve (a, ComplexMatrix (b), info); 915 return qrsolve (a, ComplexMatrix (b), info);
916 } 916 }