Changeset 2430 for trunk/sources/ALGLIB/qr.cs
- Timestamp:
- 10/15/09 13:22:07 (15 years ago)
- Location:
- trunk/sources/ALGLIB
- Files:
-
- 1 added
- 1 copied
-
. (added)
-
qr.cs (copied) (copied from trunk/sources/HeuristicLab.LinearRegression/3.2/alglib/qr.cs) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/ALGLIB/qr.cs
r2426 r2430 8 8 See subroutines comments for additional copyrights. 9 9 10 Redistribution and use in source and binary forms, with or without 11 modification, are permitted provided that the following conditions are 12 met: 13 14 - Redistributions of source code must retain the above copyright 15 notice, this list of conditions and the following disclaimer. 16 17 - Redistributions in binary form must reproduce the above copyright 18 notice, this list of conditions and the following disclaimer listed 19 in this license in the documentation and/or other materials 20 provided with the distribution. 21 22 - Neither the name of the copyright holders nor the names of its 23 contributors may be used to endorse or promote products derived from 24 this software without specific prior written permission. 25 26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 27 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 28 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 29 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 30 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 31 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 32 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 33 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 34 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 35 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 36 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 10 >>> SOURCE LICENSE >>> 11 This program is free software; you can redistribute it and/or modify 12 it under the terms of the GNU General Public License as published by 13 the Free Software Foundation (www.fsf.org); either version 2 of the 14 License, or (at your option) any later version. 15 16 This program is distributed in the hope that it will be useful, 17 but WITHOUT ANY WARRANTY; without even the implied warranty of 18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 GNU General Public License for more details. 20 21 A copy of the GNU General Public License is available at 22 http://www.fsf.org/licensing/licenses 23 24 >>> END OF LICENSE >>> 37 25 *************************************************************************/ 38 26 39 27 using System; 40 28 41 class qr 29 namespace alglib 42 30 { 43 /************************************************************************* 44 QR decomposition of a rectangular matrix of size MxN 45 46 Input parameters: 47 A - matrix A whose indexes range within [0..M-1, 0..N-1]. 48 M - number of rows in matrix A. 49 N - number of columns in matrix A. 50 51 Output parameters: 52 A - matrices Q and R in compact form (see below). 53 Tau - array of scalar factors which are used to form 54 matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)]. 55 56 Matrix A is represented as A = QR, where Q is an orthogonal matrix of size 57 MxM, R - upper triangular (or upper trapezoid) matrix of size M x N. 58 59 The elements of matrix R are located on and above the main diagonal of 60 matrix A. The elements which are located in Tau array and below the main 61 diagonal of matrix A are used to form matrix Q as follows: 62 63 Matrix Q is represented as a product of elementary reflections 64 65 Q = H(0)*H(2)*...*H(k-1), 66 67 where k = min(m,n), and each H(i) is in the form 68 69 H(i) = 1 - tau * v * (v^T) 70 71 where tau is a scalar stored in Tau[I]; v - real vector, 72 so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i). 73 74 -- LAPACK routine (version 3.0) -- 75 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 76 Courant Institute, Argonne National Lab, and Rice University 77 February 29, 1992. 78 Translation from FORTRAN to pseudocode (AlgoPascal) 79 by Sergey Bochkanov, ALGLIB project, 2005-2007. 80 *************************************************************************/ 81 public static void rmatrixqr(ref double[,] a, 82 int m, 83 int n, 84 ref double[] tau) 31 public class qr 85 32 { 86 double[] work = new double[0]; 87 double[] t = new double[0]; 88 int i = 0; 89 int k = 0; 90 int minmn = 0; 91 double tmp = 0; 92 int i_ = 0; 93 int i1_ = 0; 94 95 if( m<=0 | n<=0 ) 96 { 97 return; 98 } 99 minmn = Math.Min(m, n); 100 work = new double[n-1+1]; 101 t = new double[m+1]; 102 tau = new double[minmn-1+1]; 103 104 // 105 // Test the input arguments 106 // 107 k = minmn; 108 for(i=0; i<=k-1; i++) 109 { 110 111 // 112 // Generate elementary reflector H(i) to annihilate A(i+1:m,i) 113 // 114 i1_ = (i) - (1); 115 for(i_=1; i_<=m-i;i_++) 116 { 117 t[i_] = a[i_+i1_,i]; 118 } 119 reflections.generatereflection(ref t, m-i, ref tmp); 120 tau[i] = tmp; 121 i1_ = (1) - (i); 122 for(i_=i; i_<=m-1;i_++) 123 { 124 a[i_,i] = t[i_+i1_]; 125 } 126 t[1] = 1; 127 if( i<n ) 33 /************************************************************************* 34 QR decomposition of a rectangular matrix of size MxN 35 36 Input parameters: 37 A - matrix A whose indexes range within [0..M-1, 0..N-1]. 38 M - number of rows in matrix A. 39 N - number of columns in matrix A. 40 41 Output parameters: 42 A - matrices Q and R in compact form (see below). 43 Tau - array of scalar factors which are used to form 44 matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)]. 45 46 Matrix A is represented as A = QR, where Q is an orthogonal matrix of size 47 MxM, R - upper triangular (or upper trapezoid) matrix of size M x N. 48 49 The elements of matrix R are located on and above the main diagonal of 50 matrix A. The elements which are located in Tau array and below the main 51 diagonal of matrix A are used to form matrix Q as follows: 52 53 Matrix Q is represented as a product of elementary reflections 54 55 Q = H(0)*H(2)*...*H(k-1), 56 57 where k = min(m,n), and each H(i) is in the form 58 59 H(i) = 1 - tau * v * (v^T) 60 61 where tau is a scalar stored in Tau[I]; v - real vector, 62 so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i). 63 64 -- LAPACK routine (version 3.0) -- 65 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 66 Courant Institute, Argonne National Lab, and Rice University 67 February 29, 1992. 68 Translation from FORTRAN to pseudocode (AlgoPascal) 69 by Sergey Bochkanov, ALGLIB project, 2005-2007. 70 *************************************************************************/ 71 public static void rmatrixqr(ref double[,] a, 72 int m, 73 int n, 74 ref double[] tau) 75 { 76 double[] work = new double[0]; 77 double[] t = new double[0]; 78 int i = 0; 79 int k = 0; 80 int minmn = 0; 81 double tmp = 0; 82 int i_ = 0; 83 int i1_ = 0; 84 85 if( m<=0 | n<=0 ) 86 { 87 return; 88 } 89 minmn = Math.Min(m, n); 90 work = new double[n-1+1]; 91 t = new double[m+1]; 92 tau = new double[minmn-1+1]; 93 94 // 95 // Test the input arguments 96 // 97 k = minmn; 98 for(i=0; i<=k-1; i++) 128 99 { 129 100 130 101 // 131 // Apply H(i) to A(i:m-1,i+1:n-1) from the left 132 // 133 reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work); 134 } 135 } 136 } 137 138 139 /************************************************************************* 140 Partial unpacking of matrix Q from the QR decomposition of a matrix A 141 142 Input parameters: 143 A - matrices Q and R in compact form. 144 Output of RMatrixQR subroutine. 145 M - number of rows in given matrix A. M>=0. 146 N - number of columns in given matrix A. N>=0. 147 Tau - scalar factors which are used to form Q. 148 Output of the RMatrixQR subroutine. 149 QColumns - required number of columns of matrix Q. M>=QColumns>=0. 150 151 Output parameters: 152 Q - first QColumns columns of matrix Q. 153 Array whose indexes range within [0..M-1, 0..QColumns-1]. 154 If QColumns=0, the array remains unchanged. 155 156 -- ALGLIB -- 157 Copyright 2005 by Bochkanov Sergey 158 *************************************************************************/ 159 public static void rmatrixqrunpackq(ref double[,] a, 160 int m, 161 int n, 162 ref double[] tau, 163 int qcolumns, 164 ref double[,] q) 165 { 166 int i = 0; 167 int j = 0; 168 int k = 0; 169 int minmn = 0; 170 double[] v = new double[0]; 171 double[] work = new double[0]; 172 int i_ = 0; 173 int i1_ = 0; 174 175 System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!"); 176 if( m<=0 | n<=0 | qcolumns<=0 ) 177 { 178 return; 179 } 180 181 // 182 // init 183 // 184 minmn = Math.Min(m, n); 185 k = Math.Min(minmn, qcolumns); 186 q = new double[m-1+1, qcolumns-1+1]; 187 v = new double[m+1]; 188 work = new double[qcolumns-1+1]; 189 for(i=0; i<=m-1; i++) 190 { 191 for(j=0; j<=qcolumns-1; j++) 192 { 193 if( i==j ) 194 { 195 q[i,j] = 1; 196 } 197 else 198 { 199 q[i,j] = 0; 200 } 201 } 202 } 203 204 // 205 // unpack Q 206 // 207 for(i=k-1; i>=0; i--) 208 { 209 210 // 211 // Apply H(i) 212 // 213 i1_ = (i) - (1); 214 for(i_=1; i_<=m-i;i_++) 215 { 216 v[i_] = a[i_+i1_,i]; 217 } 218 v[1] = 1; 219 reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m-1, 0, qcolumns-1, ref work); 220 } 221 } 222 223 224 /************************************************************************* 225 Unpacking of matrix R from the QR decomposition of a matrix A 226 227 Input parameters: 228 A - matrices Q and R in compact form. 229 Output of RMatrixQR subroutine. 230 M - number of rows in given matrix A. M>=0. 231 N - number of columns in given matrix A. N>=0. 232 233 Output parameters: 234 R - matrix R, array[0..M-1, 0..N-1]. 235 236 -- ALGLIB -- 237 Copyright 2005 by Bochkanov Sergey 238 *************************************************************************/ 239 public static void rmatrixqrunpackr(ref double[,] a, 240 int m, 241 int n, 242 ref double[,] r) 243 { 244 int i = 0; 245 int k = 0; 246 int i_ = 0; 247 248 if( m<=0 | n<=0 ) 249 { 250 return; 251 } 252 k = Math.Min(m, n); 253 r = new double[m-1+1, n-1+1]; 254 for(i=0; i<=n-1; i++) 255 { 256 r[0,i] = 0; 257 } 258 for(i=1; i<=m-1; i++) 259 { 260 for(i_=0; i_<=n-1;i_++) 261 { 262 r[i,i_] = r[0,i_]; 263 } 264 } 265 for(i=0; i<=k-1; i++) 266 { 267 for(i_=i; i_<=n-1;i_++) 268 { 269 r[i,i_] = a[i,i_]; 270 } 271 } 272 } 273 274 275 /************************************************************************* 276 Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement. 277 *************************************************************************/ 278 public static void qrdecomposition(ref double[,] a, 279 int m, 280 int n, 281 ref double[] tau) 282 { 283 double[] work = new double[0]; 284 double[] t = new double[0]; 285 int i = 0; 286 int k = 0; 287 int mmip1 = 0; 288 int minmn = 0; 289 double tmp = 0; 290 int i_ = 0; 291 int i1_ = 0; 292 293 minmn = Math.Min(m, n); 294 work = new double[n+1]; 295 t = new double[m+1]; 296 tau = new double[minmn+1]; 297 298 // 299 // Test the input arguments 300 // 301 k = Math.Min(m, n); 302 for(i=1; i<=k; i++) 303 { 304 305 // 306 // Generate elementary reflector H(i) to annihilate A(i+1:m,i) 307 // 308 mmip1 = m-i+1; 309 i1_ = (i) - (1); 310 for(i_=1; i_<=mmip1;i_++) 311 { 312 t[i_] = a[i_+i1_,i]; 313 } 314 reflections.generatereflection(ref t, mmip1, ref tmp); 315 tau[i] = tmp; 316 i1_ = (1) - (i); 317 for(i_=i; i_<=m;i_++) 318 { 319 a[i_,i] = t[i_+i1_]; 320 } 321 t[1] = 1; 322 if( i<n ) 102 // Generate elementary reflector H(i) to annihilate A(i+1:m,i) 103 // 104 i1_ = (i) - (1); 105 for(i_=1; i_<=m-i;i_++) 106 { 107 t[i_] = a[i_+i1_,i]; 108 } 109 reflections.generatereflection(ref t, m-i, ref tmp); 110 tau[i] = tmp; 111 i1_ = (1) - (i); 112 for(i_=i; i_<=m-1;i_++) 113 { 114 a[i_,i] = t[i_+i1_]; 115 } 116 t[1] = 1; 117 if( i<n ) 118 { 119 120 // 121 // Apply H(i) to A(i:m-1,i+1:n-1) from the left 122 // 123 reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work); 124 } 125 } 126 } 127 128 129 /************************************************************************* 130 Partial unpacking of matrix Q from the QR decomposition of a matrix A 131 132 Input parameters: 133 A - matrices Q and R in compact form. 134 Output of RMatrixQR subroutine. 135 M - number of rows in given matrix A. M>=0. 136 N - number of columns in given matrix A. N>=0. 137 Tau - scalar factors which are used to form Q. 138 Output of the RMatrixQR subroutine. 139 QColumns - required number of columns of matrix Q. M>=QColumns>=0. 140 141 Output parameters: 142 Q - first QColumns columns of matrix Q. 143 Array whose indexes range within [0..M-1, 0..QColumns-1]. 144 If QColumns=0, the array remains unchanged. 145 146 -- ALGLIB -- 147 Copyright 2005 by Bochkanov Sergey 148 *************************************************************************/ 149 public static void rmatrixqrunpackq(ref double[,] a, 150 int m, 151 int n, 152 ref double[] tau, 153 int qcolumns, 154 ref double[,] q) 155 { 156 int i = 0; 157 int j = 0; 158 int k = 0; 159 int minmn = 0; 160 double[] v = new double[0]; 161 double[] work = new double[0]; 162 int i_ = 0; 163 int i1_ = 0; 164 165 System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!"); 166 if( m<=0 | n<=0 | qcolumns<=0 ) 167 { 168 return; 169 } 170 171 // 172 // init 173 // 174 minmn = Math.Min(m, n); 175 k = Math.Min(minmn, qcolumns); 176 q = new double[m-1+1, qcolumns-1+1]; 177 v = new double[m+1]; 178 work = new double[qcolumns-1+1]; 179 for(i=0; i<=m-1; i++) 180 { 181 for(j=0; j<=qcolumns-1; j++) 182 { 183 if( i==j ) 184 { 185 q[i,j] = 1; 186 } 187 else 188 { 189 q[i,j] = 0; 190 } 191 } 192 } 193 194 // 195 // unpack Q 196 // 197 for(i=k-1; i>=0; i--) 323 198 { 324 199 325 200 // 326 // Apply H(i) to A(i:m,i+1:n) from the left 327 // 328 reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work); 329 } 330 } 331 } 332 333 334 /************************************************************************* 335 Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement. 336 *************************************************************************/ 337 public static void unpackqfromqr(ref double[,] a, 338 int m, 339 int n, 340 ref double[] tau, 341 int qcolumns, 342 ref double[,] q) 343 { 344 int i = 0; 345 int j = 0; 346 int k = 0; 347 int minmn = 0; 348 double[] v = new double[0]; 349 double[] work = new double[0]; 350 int vm = 0; 351 int i_ = 0; 352 int i1_ = 0; 353 354 System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!"); 355 if( m==0 | n==0 | qcolumns==0 ) 356 { 357 return; 358 } 359 360 // 361 // init 362 // 363 minmn = Math.Min(m, n); 364 k = Math.Min(minmn, qcolumns); 365 q = new double[m+1, qcolumns+1]; 366 v = new double[m+1]; 367 work = new double[qcolumns+1]; 368 for(i=1; i<=m; i++) 369 { 370 for(j=1; j<=qcolumns; j++) 371 { 372 if( i==j ) 373 { 374 q[i,j] = 1; 375 } 376 else 377 { 378 q[i,j] = 0; 379 } 380 } 381 } 382 383 // 384 // unpack Q 385 // 386 for(i=k; i>=1; i--) 387 { 388 389 // 390 // Apply H(i) 391 // 392 vm = m-i+1; 393 i1_ = (i) - (1); 394 for(i_=1; i_<=vm;i_++) 395 { 396 v[i_] = a[i_+i1_,i]; 397 } 398 v[1] = 1; 399 reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work); 400 } 401 } 402 403 404 /************************************************************************* 405 Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement. 406 *************************************************************************/ 407 public static void qrdecompositionunpacked(double[,] a, 408 int m, 409 int n, 410 ref double[,] q, 411 ref double[,] r) 412 { 413 int i = 0; 414 int k = 0; 415 double[] tau = new double[0]; 416 double[] work = new double[0]; 417 double[] v = new double[0]; 418 int i_ = 0; 419 420 a = (double[,])a.Clone(); 421 422 k = Math.Min(m, n); 423 if( n<=0 ) 424 { 425 return; 426 } 427 work = new double[m+1]; 428 v = new double[m+1]; 429 q = new double[m+1, m+1]; 430 r = new double[m+1, n+1]; 431 432 // 433 // QRDecomposition 434 // 435 qrdecomposition(ref a, m, n, ref tau); 436 437 // 438 // R 439 // 440 for(i=1; i<=n; i++) 441 { 442 r[1,i] = 0; 443 } 444 for(i=2; i<=m; i++) 445 { 446 for(i_=1; i_<=n;i_++) 447 { 448 r[i,i_] = r[1,i_]; 449 } 450 } 451 for(i=1; i<=k; i++) 452 { 453 for(i_=i; i_<=n;i_++) 454 { 455 r[i,i_] = a[i,i_]; 456 } 457 } 458 459 // 460 // Q 461 // 462 unpackqfromqr(ref a, m, n, ref tau, m, ref q); 201 // Apply H(i) 202 // 203 i1_ = (i) - (1); 204 for(i_=1; i_<=m-i;i_++) 205 { 206 v[i_] = a[i_+i1_,i]; 207 } 208 v[1] = 1; 209 reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m-1, 0, qcolumns-1, ref work); 210 } 211 } 212 213 214 /************************************************************************* 215 Unpacking of matrix R from the QR decomposition of a matrix A 216 217 Input parameters: 218 A - matrices Q and R in compact form. 219 Output of RMatrixQR subroutine. 220 M - number of rows in given matrix A. M>=0. 221 N - number of columns in given matrix A. N>=0. 222 223 Output parameters: 224 R - matrix R, array[0..M-1, 0..N-1]. 225 226 -- ALGLIB -- 227 Copyright 2005 by Bochkanov Sergey 228 *************************************************************************/ 229 public static void rmatrixqrunpackr(ref double[,] a, 230 int m, 231 int n, 232 ref double[,] r) 233 { 234 int i = 0; 235 int k = 0; 236 int i_ = 0; 237 238 if( m<=0 | n<=0 ) 239 { 240 return; 241 } 242 k = Math.Min(m, n); 243 r = new double[m-1+1, n-1+1]; 244 for(i=0; i<=n-1; i++) 245 { 246 r[0,i] = 0; 247 } 248 for(i=1; i<=m-1; i++) 249 { 250 for(i_=0; i_<=n-1;i_++) 251 { 252 r[i,i_] = r[0,i_]; 253 } 254 } 255 for(i=0; i<=k-1; i++) 256 { 257 for(i_=i; i_<=n-1;i_++) 258 { 259 r[i,i_] = a[i,i_]; 260 } 261 } 262 } 263 264 265 /************************************************************************* 266 Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement. 267 *************************************************************************/ 268 public static void qrdecomposition(ref double[,] a, 269 int m, 270 int n, 271 ref double[] tau) 272 { 273 double[] work = new double[0]; 274 double[] t = new double[0]; 275 int i = 0; 276 int k = 0; 277 int mmip1 = 0; 278 int minmn = 0; 279 double tmp = 0; 280 int i_ = 0; 281 int i1_ = 0; 282 283 minmn = Math.Min(m, n); 284 work = new double[n+1]; 285 t = new double[m+1]; 286 tau = new double[minmn+1]; 287 288 // 289 // Test the input arguments 290 // 291 k = Math.Min(m, n); 292 for(i=1; i<=k; i++) 293 { 294 295 // 296 // Generate elementary reflector H(i) to annihilate A(i+1:m,i) 297 // 298 mmip1 = m-i+1; 299 i1_ = (i) - (1); 300 for(i_=1; i_<=mmip1;i_++) 301 { 302 t[i_] = a[i_+i1_,i]; 303 } 304 reflections.generatereflection(ref t, mmip1, ref tmp); 305 tau[i] = tmp; 306 i1_ = (1) - (i); 307 for(i_=i; i_<=m;i_++) 308 { 309 a[i_,i] = t[i_+i1_]; 310 } 311 t[1] = 1; 312 if( i<n ) 313 { 314 315 // 316 // Apply H(i) to A(i:m,i+1:n) from the left 317 // 318 reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work); 319 } 320 } 321 } 322 323 324 /************************************************************************* 325 Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement. 326 *************************************************************************/ 327 public static void unpackqfromqr(ref double[,] a, 328 int m, 329 int n, 330 ref double[] tau, 331 int qcolumns, 332 ref double[,] q) 333 { 334 int i = 0; 335 int j = 0; 336 int k = 0; 337 int minmn = 0; 338 double[] v = new double[0]; 339 double[] work = new double[0]; 340 int vm = 0; 341 int i_ = 0; 342 int i1_ = 0; 343 344 System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!"); 345 if( m==0 | n==0 | qcolumns==0 ) 346 { 347 return; 348 } 349 350 // 351 // init 352 // 353 minmn = Math.Min(m, n); 354 k = Math.Min(minmn, qcolumns); 355 q = new double[m+1, qcolumns+1]; 356 v = new double[m+1]; 357 work = new double[qcolumns+1]; 358 for(i=1; i<=m; i++) 359 { 360 for(j=1; j<=qcolumns; j++) 361 { 362 if( i==j ) 363 { 364 q[i,j] = 1; 365 } 366 else 367 { 368 q[i,j] = 0; 369 } 370 } 371 } 372 373 // 374 // unpack Q 375 // 376 for(i=k; i>=1; i--) 377 { 378 379 // 380 // Apply H(i) 381 // 382 vm = m-i+1; 383 i1_ = (i) - (1); 384 for(i_=1; i_<=vm;i_++) 385 { 386 v[i_] = a[i_+i1_,i]; 387 } 388 v[1] = 1; 389 reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work); 390 } 391 } 392 393 394 /************************************************************************* 395 Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement. 396 *************************************************************************/ 397 public static void qrdecompositionunpacked(double[,] a, 398 int m, 399 int n, 400 ref double[,] q, 401 ref double[,] r) 402 { 403 int i = 0; 404 int k = 0; 405 double[] tau = new double[0]; 406 double[] work = new double[0]; 407 double[] v = new double[0]; 408 int i_ = 0; 409 410 a = (double[,])a.Clone(); 411 412 k = Math.Min(m, n); 413 if( n<=0 ) 414 { 415 return; 416 } 417 work = new double[m+1]; 418 v = new double[m+1]; 419 q = new double[m+1, m+1]; 420 r = new double[m+1, n+1]; 421 422 // 423 // QRDecomposition 424 // 425 qrdecomposition(ref a, m, n, ref tau); 426 427 // 428 // R 429 // 430 for(i=1; i<=n; i++) 431 { 432 r[1,i] = 0; 433 } 434 for(i=2; i<=m; i++) 435 { 436 for(i_=1; i_<=n;i_++) 437 { 438 r[i,i_] = r[1,i_]; 439 } 440 } 441 for(i=1; i<=k; i++) 442 { 443 for(i_=i; i_<=n;i_++) 444 { 445 r[i,i_] = a[i,i_]; 446 } 447 } 448 449 // 450 // Q 451 // 452 unpackqfromqr(ref a, m, n, ref tau, m, ref q); 453 } 463 454 } 464 455 }
Note: See TracChangeset
for help on using the changeset viewer.
