Changeset 4068 for trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/spdgevd.cs
- Timestamp:
- 07/22/10 00:44:01 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/spdgevd.cs
r3839 r4068 19 19 *************************************************************************/ 20 20 21 using System; 22 23 namespace alglib 24 { 25 public class spdgevd 26 { 27 /************************************************************************* 28 Algorithm for solving the following generalized symmetric positive-definite 29 eigenproblem: 30 A*x = lambda*B*x (1) or 31 A*B*x = lambda*x (2) or 32 B*A*x = lambda*x (3). 33 where A is a symmetric matrix, B - symmetric positive-definite matrix. 34 The problem is solved by reducing it to an ordinary symmetric eigenvalue 35 problem. 36 37 Input parameters: 38 A - symmetric matrix which is given by its upper or lower 39 triangular part. 40 Array whose indexes range within [0..N-1, 0..N-1]. 41 N - size of matrices A and B. 42 IsUpperA - storage format of matrix A. 43 B - symmetric positive-definite matrix which is given by 44 its upper or lower triangular part. 45 Array whose indexes range within [0..N-1, 0..N-1]. 46 IsUpperB - storage format of matrix B. 47 ZNeeded - if ZNeeded is equal to: 48 * 0, the eigenvectors are not returned; 49 * 1, the eigenvectors are returned. 50 ProblemType - if ProblemType is equal to: 51 * 1, the following problem is solved: A*x = lambda*B*x; 52 * 2, the following problem is solved: A*B*x = lambda*x; 53 * 3, the following problem is solved: B*A*x = lambda*x. 54 55 Output parameters: 56 D - eigenvalues in ascending order. 57 Array whose index ranges within [0..N-1]. 58 Z - if ZNeeded is equal to: 59 * 0, Z hasnt changed; 60 * 1, Z contains eigenvectors. 61 Array whose indexes range within [0..N-1, 0..N-1]. 62 The eigenvectors are stored in matrix columns. It should 63 be noted that the eigenvectors in such problems do not 64 form an orthogonal system. 65 66 Result: 67 True, if the problem was solved successfully. 68 False, if the error occurred during the Cholesky decomposition of matrix 69 B (the matrix isnt positive-definite) or during the work of the iterative 70 algorithm for solving the symmetric eigenproblem. 71 72 See also the GeneralizedSymmetricDefiniteEVDReduce subroutine. 73 74 -- ALGLIB -- 75 Copyright 1.28.2006 by Bochkanov Sergey 76 *************************************************************************/ 77 public static bool smatrixgevd(double[,] a, 78 int n, 79 bool isuppera, 80 ref double[,] b, 81 bool isupperb, 82 int zneeded, 83 int problemtype, 84 ref double[] d, 85 ref double[,] z) 86 { 87 bool result = new bool(); 88 double[,] r = new double[0,0]; 89 double[,] t = new double[0,0]; 90 bool isupperr = new bool(); 91 int j1 = 0; 92 int j2 = 0; 93 int j1inc = 0; 94 int j2inc = 0; 95 int i = 0; 96 int j = 0; 97 double v = 0; 98 int i_ = 0; 99 100 a = (double[,])a.Clone(); 101 102 103 // 104 // Reduce and solve 105 // 106 result = smatrixgevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr); 107 if( !result ) 108 { 109 return result; 110 } 111 result = evd.smatrixevd(a, n, zneeded, isuppera, ref d, ref t); 112 if( !result ) 113 { 114 return result; 115 } 116 117 // 118 // Transform eigenvectors if needed 119 // 120 if( zneeded!=0 ) 121 { 122 123 // 124 // fill Z with zeros 125 // 126 z = new double[n-1+1, n-1+1]; 127 for(j=0; j<=n-1; j++) 128 { 129 z[0,j] = 0.0; 130 } 131 for(i=1; i<=n-1; i++) 132 { 133 for(i_=0; i_<=n-1;i_++) 134 { 135 z[i,i_] = z[0,i_]; 136 } 137 } 138 139 // 140 // Setup R properties 141 // 142 if( isupperr ) 143 { 144 j1 = 0; 145 j2 = n-1; 146 j1inc = +1; 147 j2inc = 0; 148 } 149 else 150 { 151 j1 = 0; 152 j2 = 0; 153 j1inc = 0; 154 j2inc = +1; 155 } 156 157 // 158 // Calculate R*Z 159 // 160 for(i=0; i<=n-1; i++) 161 { 162 for(j=j1; j<=j2; j++) 163 { 164 v = r[i,j]; 165 for(i_=0; i_<=n-1;i_++) 166 { 167 z[i,i_] = z[i,i_] + v*t[j,i_]; 168 } 169 } 170 j1 = j1+j1inc; 171 j2 = j2+j2inc; 172 } 173 } 21 22 namespace alglib { 23 public class spdgevd { 24 /************************************************************************* 25 Algorithm for solving the following generalized symmetric positive-definite 26 eigenproblem: 27 A*x = lambda*B*x (1) or 28 A*B*x = lambda*x (2) or 29 B*A*x = lambda*x (3). 30 where A is a symmetric matrix, B - symmetric positive-definite matrix. 31 The problem is solved by reducing it to an ordinary symmetric eigenvalue 32 problem. 33 34 Input parameters: 35 A - symmetric matrix which is given by its upper or lower 36 triangular part. 37 Array whose indexes range within [0..N-1, 0..N-1]. 38 N - size of matrices A and B. 39 IsUpperA - storage format of matrix A. 40 B - symmetric positive-definite matrix which is given by 41 its upper or lower triangular part. 42 Array whose indexes range within [0..N-1, 0..N-1]. 43 IsUpperB - storage format of matrix B. 44 ZNeeded - if ZNeeded is equal to: 45 * 0, the eigenvectors are not returned; 46 * 1, the eigenvectors are returned. 47 ProblemType - if ProblemType is equal to: 48 * 1, the following problem is solved: A*x = lambda*B*x; 49 * 2, the following problem is solved: A*B*x = lambda*x; 50 * 3, the following problem is solved: B*A*x = lambda*x. 51 52 Output parameters: 53 D - eigenvalues in ascending order. 54 Array whose index ranges within [0..N-1]. 55 Z - if ZNeeded is equal to: 56 * 0, Z hasnt changed; 57 * 1, Z contains eigenvectors. 58 Array whose indexes range within [0..N-1, 0..N-1]. 59 The eigenvectors are stored in matrix columns. It should 60 be noted that the eigenvectors in such problems do not 61 form an orthogonal system. 62 63 Result: 64 True, if the problem was solved successfully. 65 False, if the error occurred during the Cholesky decomposition of matrix 66 B (the matrix isnt positive-definite) or during the work of the iterative 67 algorithm for solving the symmetric eigenproblem. 68 69 See also the GeneralizedSymmetricDefiniteEVDReduce subroutine. 70 71 -- ALGLIB -- 72 Copyright 1.28.2006 by Bochkanov Sergey 73 *************************************************************************/ 74 public static bool smatrixgevd(double[,] a, 75 int n, 76 bool isuppera, 77 ref double[,] b, 78 bool isupperb, 79 int zneeded, 80 int problemtype, 81 ref double[] d, 82 ref double[,] z) { 83 bool result = new bool(); 84 double[,] r = new double[0, 0]; 85 double[,] t = new double[0, 0]; 86 bool isupperr = new bool(); 87 int j1 = 0; 88 int j2 = 0; 89 int j1inc = 0; 90 int j2inc = 0; 91 int i = 0; 92 int j = 0; 93 double v = 0; 94 int i_ = 0; 95 96 a = (double[,])a.Clone(); 97 98 99 // 100 // Reduce and solve 101 // 102 result = smatrixgevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr); 103 if (!result) { 104 return result; 105 } 106 result = evd.smatrixevd(a, n, zneeded, isuppera, ref d, ref t); 107 if (!result) { 108 return result; 109 } 110 111 // 112 // Transform eigenvectors if needed 113 // 114 if (zneeded != 0) { 115 116 // 117 // fill Z with zeros 118 // 119 z = new double[n - 1 + 1, n - 1 + 1]; 120 for (j = 0; j <= n - 1; j++) { 121 z[0, j] = 0.0; 122 } 123 for (i = 1; i <= n - 1; i++) { 124 for (i_ = 0; i_ <= n - 1; i_++) { 125 z[i, i_] = z[0, i_]; 126 } 127 } 128 129 // 130 // Setup R properties 131 // 132 if (isupperr) { 133 j1 = 0; 134 j2 = n - 1; 135 j1inc = +1; 136 j2inc = 0; 137 } else { 138 j1 = 0; 139 j2 = 0; 140 j1inc = 0; 141 j2inc = +1; 142 } 143 144 // 145 // Calculate R*Z 146 // 147 for (i = 0; i <= n - 1; i++) { 148 for (j = j1; j <= j2; j++) { 149 v = r[i, j]; 150 for (i_ = 0; i_ <= n - 1; i_++) { 151 z[i, i_] = z[i, i_] + v * t[j, i_]; 152 } 153 } 154 j1 = j1 + j1inc; 155 j2 = j2 + j2inc; 156 } 157 } 158 return result; 159 } 160 161 162 /************************************************************************* 163 Algorithm for reduction of the following generalized symmetric positive- 164 definite eigenvalue problem: 165 A*x = lambda*B*x (1) or 166 A*B*x = lambda*x (2) or 167 B*A*x = lambda*x (3) 168 to the symmetric eigenvalues problem C*y = lambda*y (eigenvalues of this and 169 the given problems are the same, and the eigenvectors of the given problem 170 could be obtained by multiplying the obtained eigenvectors by the 171 transformation matrix x = R*y). 172 173 Here A is a symmetric matrix, B - symmetric positive-definite matrix. 174 175 Input parameters: 176 A - symmetric matrix which is given by its upper or lower 177 triangular part. 178 Array whose indexes range within [0..N-1, 0..N-1]. 179 N - size of matrices A and B. 180 IsUpperA - storage format of matrix A. 181 B - symmetric positive-definite matrix which is given by 182 its upper or lower triangular part. 183 Array whose indexes range within [0..N-1, 0..N-1]. 184 IsUpperB - storage format of matrix B. 185 ProblemType - if ProblemType is equal to: 186 * 1, the following problem is solved: A*x = lambda*B*x; 187 * 2, the following problem is solved: A*B*x = lambda*x; 188 * 3, the following problem is solved: B*A*x = lambda*x. 189 190 Output parameters: 191 A - symmetric matrix which is given by its upper or lower 192 triangle depending on IsUpperA. Contains matrix C. 193 Array whose indexes range within [0..N-1, 0..N-1]. 194 R - upper triangular or low triangular transformation matrix 195 which is used to obtain the eigenvectors of a given problem 196 as the product of eigenvectors of C (from the right) and 197 matrix R (from the left). If the matrix is upper 198 triangular, the elements below the main diagonal 199 are equal to 0 (and vice versa). Thus, we can perform 200 the multiplication without taking into account the 201 internal structure (which is an easier though less 202 effective way). 203 Array whose indexes range within [0..N-1, 0..N-1]. 204 IsUpperR - type of matrix R (upper or lower triangular). 205 206 Result: 207 True, if the problem was reduced successfully. 208 False, if the error occurred during the Cholesky decomposition of 209 matrix B (the matrix is not positive-definite). 210 211 -- ALGLIB -- 212 Copyright 1.28.2006 by Bochkanov Sergey 213 *************************************************************************/ 214 public static bool smatrixgevdreduce(ref double[,] a, 215 int n, 216 bool isuppera, 217 ref double[,] b, 218 bool isupperb, 219 int problemtype, 220 ref double[,] r, 221 ref bool isupperr) { 222 bool result = new bool(); 223 double[,] t = new double[0, 0]; 224 double[] w1 = new double[0]; 225 double[] w2 = new double[0]; 226 double[] w3 = new double[0]; 227 int i = 0; 228 int j = 0; 229 double v = 0; 230 matinv.matinvreport rep = new matinv.matinvreport(); 231 int info = 0; 232 int i_ = 0; 233 int i1_ = 0; 234 235 System.Diagnostics.Debug.Assert(n > 0, "SMatrixGEVDReduce: N<=0!"); 236 System.Diagnostics.Debug.Assert(problemtype == 1 | problemtype == 2 | problemtype == 3, "SMatrixGEVDReduce: incorrect ProblemType!"); 237 result = true; 238 239 // 240 // Problem 1: A*x = lambda*B*x 241 // 242 // Reducing to: 243 // C*y = lambda*y 244 // C = L^(-1) * A * L^(-T) 245 // x = L^(-T) * y 246 // 247 if (problemtype == 1) { 248 249 // 250 // Factorize B in T: B = LL' 251 // 252 t = new double[n - 1 + 1, n - 1 + 1]; 253 if (isupperb) { 254 for (i = 0; i <= n - 1; i++) { 255 for (i_ = i; i_ <= n - 1; i_++) { 256 t[i_, i] = b[i, i_]; 257 } 258 } 259 } else { 260 for (i = 0; i <= n - 1; i++) { 261 for (i_ = 0; i_ <= i; i_++) { 262 t[i, i_] = b[i, i_]; 263 } 264 } 265 } 266 if (!trfac.spdmatrixcholesky(ref t, n, false)) { 267 result = false; 268 return result; 269 } 270 271 // 272 // Invert L in T 273 // 274 matinv.rmatrixtrinverse(ref t, n, false, false, ref info, ref rep); 275 if (info <= 0) { 276 result = false; 277 return result; 278 } 279 280 // 281 // Build L^(-1) * A * L^(-T) in R 282 // 283 w1 = new double[n + 1]; 284 w2 = new double[n + 1]; 285 r = new double[n - 1 + 1, n - 1 + 1]; 286 for (j = 1; j <= n; j++) { 287 288 // 289 // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T)) 290 // 291 i1_ = (0) - (1); 292 for (i_ = 1; i_ <= j; i_++) { 293 w1[i_] = t[j - 1, i_ + i1_]; 294 } 295 sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j - 1, ref w1, 1.0, ref w2); 296 if (isuppera) { 297 blas.matrixvectormultiply(ref a, 0, j - 1, j, n - 1, true, ref w1, 1, j, 1.0, ref w2, j + 1, n, 0.0); 298 } else { 299 blas.matrixvectormultiply(ref a, j, n - 1, 0, j - 1, false, ref w1, 1, j, 1.0, ref w2, j + 1, n, 0.0); 300 } 301 302 // 303 // Form l(i)*w2 (here l(i) is i-th row of L^(-1)) 304 // 305 for (i = 1; i <= n; i++) { 306 i1_ = (1) - (0); 307 v = 0.0; 308 for (i_ = 0; i_ <= i - 1; i_++) { 309 v += t[i - 1, i_] * w2[i_ + i1_]; 310 } 311 r[i - 1, j - 1] = v; 312 } 313 } 314 315 // 316 // Copy R to A 317 // 318 for (i = 0; i <= n - 1; i++) { 319 for (i_ = 0; i_ <= n - 1; i_++) { 320 a[i, i_] = r[i, i_]; 321 } 322 } 323 324 // 325 // Copy L^(-1) from T to R and transpose 326 // 327 isupperr = true; 328 for (i = 0; i <= n - 1; i++) { 329 for (j = 0; j <= i - 1; j++) { 330 r[i, j] = 0; 331 } 332 } 333 for (i = 0; i <= n - 1; i++) { 334 for (i_ = i; i_ <= n - 1; i_++) { 335 r[i, i_] = t[i_, i]; 336 } 337 } 338 return result; 339 } 340 341 // 342 // Problem 2: A*B*x = lambda*x 343 // or 344 // problem 3: B*A*x = lambda*x 345 // 346 // Reducing to: 347 // C*y = lambda*y 348 // C = U * A * U' 349 // B = U'* U 350 // 351 if (problemtype == 2 | problemtype == 3) { 352 353 // 354 // Factorize B in T: B = U'*U 355 // 356 t = new double[n - 1 + 1, n - 1 + 1]; 357 if (isupperb) { 358 for (i = 0; i <= n - 1; i++) { 359 for (i_ = i; i_ <= n - 1; i_++) { 360 t[i, i_] = b[i, i_]; 361 } 362 } 363 } else { 364 for (i = 0; i <= n - 1; i++) { 365 for (i_ = i; i_ <= n - 1; i_++) { 366 t[i, i_] = b[i_, i]; 367 } 368 } 369 } 370 if (!trfac.spdmatrixcholesky(ref t, n, true)) { 371 result = false; 372 return result; 373 } 374 375 // 376 // Build U * A * U' in R 377 // 378 w1 = new double[n + 1]; 379 w2 = new double[n + 1]; 380 w3 = new double[n + 1]; 381 r = new double[n - 1 + 1, n - 1 + 1]; 382 for (j = 1; j <= n; j++) { 383 384 // 385 // Form w2 = A * u'(j) (here u'(j) is j-th column of U') 386 // 387 i1_ = (j - 1) - (1); 388 for (i_ = 1; i_ <= n - j + 1; i_++) { 389 w1[i_] = t[j - 1, i_ + i1_]; 390 } 391 sblas.symmetricmatrixvectormultiply(ref a, isuppera, j - 1, n - 1, ref w1, 1.0, ref w3); 392 i1_ = (1) - (j); 393 for (i_ = j; i_ <= n; i_++) { 394 w2[i_] = w3[i_ + i1_]; 395 } 396 i1_ = (j - 1) - (j); 397 for (i_ = j; i_ <= n; i_++) { 398 w1[i_] = t[j - 1, i_ + i1_]; 399 } 400 if (isuppera) { 401 blas.matrixvectormultiply(ref a, 0, j - 2, j - 1, n - 1, false, ref w1, j, n, 1.0, ref w2, 1, j - 1, 0.0); 402 } else { 403 blas.matrixvectormultiply(ref a, j - 1, n - 1, 0, j - 2, true, ref w1, j, n, 1.0, ref w2, 1, j - 1, 0.0); 404 } 405 406 // 407 // Form u(i)*w2 (here u(i) is i-th row of U) 408 // 409 for (i = 1; i <= n; i++) { 410 i1_ = (i) - (i - 1); 411 v = 0.0; 412 for (i_ = i - 1; i_ <= n - 1; i_++) { 413 v += t[i - 1, i_] * w2[i_ + i1_]; 414 } 415 r[i - 1, j - 1] = v; 416 } 417 } 418 419 // 420 // Copy R to A 421 // 422 for (i = 0; i <= n - 1; i++) { 423 for (i_ = 0; i_ <= n - 1; i_++) { 424 a[i, i_] = r[i, i_]; 425 } 426 } 427 if (problemtype == 2) { 428 429 // 430 // Invert U in T 431 // 432 matinv.rmatrixtrinverse(ref t, n, true, false, ref info, ref rep); 433 if (info <= 0) { 434 result = false; 174 435 return result; 175 } 176 177 178 /************************************************************************* 179 Algorithm for reduction of the following generalized symmetric positive- 180 definite eigenvalue problem: 181 A*x = lambda*B*x (1) or 182 A*B*x = lambda*x (2) or 183 B*A*x = lambda*x (3) 184 to the symmetric eigenvalues problem C*y = lambda*y (eigenvalues of this and 185 the given problems are the same, and the eigenvectors of the given problem 186 could be obtained by multiplying the obtained eigenvectors by the 187 transformation matrix x = R*y). 188 189 Here A is a symmetric matrix, B - symmetric positive-definite matrix. 190 191 Input parameters: 192 A - symmetric matrix which is given by its upper or lower 193 triangular part. 194 Array whose indexes range within [0..N-1, 0..N-1]. 195 N - size of matrices A and B. 196 IsUpperA - storage format of matrix A. 197 B - symmetric positive-definite matrix which is given by 198 its upper or lower triangular part. 199 Array whose indexes range within [0..N-1, 0..N-1]. 200 IsUpperB - storage format of matrix B. 201 ProblemType - if ProblemType is equal to: 202 * 1, the following problem is solved: A*x = lambda*B*x; 203 * 2, the following problem is solved: A*B*x = lambda*x; 204 * 3, the following problem is solved: B*A*x = lambda*x. 205 206 Output parameters: 207 A - symmetric matrix which is given by its upper or lower 208 triangle depending on IsUpperA. Contains matrix C. 209 Array whose indexes range within [0..N-1, 0..N-1]. 210 R - upper triangular or low triangular transformation matrix 211 which is used to obtain the eigenvectors of a given problem 212 as the product of eigenvectors of C (from the right) and 213 matrix R (from the left). If the matrix is upper 214 triangular, the elements below the main diagonal 215 are equal to 0 (and vice versa). Thus, we can perform 216 the multiplication without taking into account the 217 internal structure (which is an easier though less 218 effective way). 219 Array whose indexes range within [0..N-1, 0..N-1]. 220 IsUpperR - type of matrix R (upper or lower triangular). 221 222 Result: 223 True, if the problem was reduced successfully. 224 False, if the error occurred during the Cholesky decomposition of 225 matrix B (the matrix is not positive-definite). 226 227 -- ALGLIB -- 228 Copyright 1.28.2006 by Bochkanov Sergey 229 *************************************************************************/ 230 public static bool smatrixgevdreduce(ref double[,] a, 231 int n, 232 bool isuppera, 233 ref double[,] b, 234 bool isupperb, 235 int problemtype, 236 ref double[,] r, 237 ref bool isupperr) 238 { 239 bool result = new bool(); 240 double[,] t = new double[0,0]; 241 double[] w1 = new double[0]; 242 double[] w2 = new double[0]; 243 double[] w3 = new double[0]; 244 int i = 0; 245 int j = 0; 246 double v = 0; 247 matinv.matinvreport rep = new matinv.matinvreport(); 248 int info = 0; 249 int i_ = 0; 250 int i1_ = 0; 251 252 System.Diagnostics.Debug.Assert(n>0, "SMatrixGEVDReduce: N<=0!"); 253 System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "SMatrixGEVDReduce: incorrect ProblemType!"); 254 result = true; 255 256 // 257 // Problem 1: A*x = lambda*B*x 258 // 259 // Reducing to: 260 // C*y = lambda*y 261 // C = L^(-1) * A * L^(-T) 262 // x = L^(-T) * y 263 // 264 if( problemtype==1 ) 265 { 266 267 // 268 // Factorize B in T: B = LL' 269 // 270 t = new double[n-1+1, n-1+1]; 271 if( isupperb ) 272 { 273 for(i=0; i<=n-1; i++) 274 { 275 for(i_=i; i_<=n-1;i_++) 276 { 277 t[i_,i] = b[i,i_]; 278 } 279 } 280 } 281 else 282 { 283 for(i=0; i<=n-1; i++) 284 { 285 for(i_=0; i_<=i;i_++) 286 { 287 t[i,i_] = b[i,i_]; 288 } 289 } 290 } 291 if( !trfac.spdmatrixcholesky(ref t, n, false) ) 292 { 293 result = false; 294 return result; 295 } 296 297 // 298 // Invert L in T 299 // 300 matinv.rmatrixtrinverse(ref t, n, false, false, ref info, ref rep); 301 if( info<=0 ) 302 { 303 result = false; 304 return result; 305 } 306 307 // 308 // Build L^(-1) * A * L^(-T) in R 309 // 310 w1 = new double[n+1]; 311 w2 = new double[n+1]; 312 r = new double[n-1+1, n-1+1]; 313 for(j=1; j<=n; j++) 314 { 315 316 // 317 // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T)) 318 // 319 i1_ = (0) - (1); 320 for(i_=1; i_<=j;i_++) 321 { 322 w1[i_] = t[j-1,i_+i1_]; 323 } 324 sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j-1, ref w1, 1.0, ref w2); 325 if( isuppera ) 326 { 327 blas.matrixvectormultiply(ref a, 0, j-1, j, n-1, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0); 328 } 329 else 330 { 331 blas.matrixvectormultiply(ref a, j, n-1, 0, j-1, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0); 332 } 333 334 // 335 // Form l(i)*w2 (here l(i) is i-th row of L^(-1)) 336 // 337 for(i=1; i<=n; i++) 338 { 339 i1_ = (1)-(0); 340 v = 0.0; 341 for(i_=0; i_<=i-1;i_++) 342 { 343 v += t[i-1,i_]*w2[i_+i1_]; 344 } 345 r[i-1,j-1] = v; 346 } 347 } 348 349 // 350 // Copy R to A 351 // 352 for(i=0; i<=n-1; i++) 353 { 354 for(i_=0; i_<=n-1;i_++) 355 { 356 a[i,i_] = r[i,i_]; 357 } 358 } 359 360 // 361 // Copy L^(-1) from T to R and transpose 362 // 363 isupperr = true; 364 for(i=0; i<=n-1; i++) 365 { 366 for(j=0; j<=i-1; j++) 367 { 368 r[i,j] = 0; 369 } 370 } 371 for(i=0; i<=n-1; i++) 372 { 373 for(i_=i; i_<=n-1;i_++) 374 { 375 r[i,i_] = t[i_,i]; 376 } 377 } 378 return result; 379 } 380 381 // 382 // Problem 2: A*B*x = lambda*x 383 // or 384 // problem 3: B*A*x = lambda*x 385 // 386 // Reducing to: 387 // C*y = lambda*y 388 // C = U * A * U' 389 // B = U'* U 390 // 391 if( problemtype==2 | problemtype==3 ) 392 { 393 394 // 395 // Factorize B in T: B = U'*U 396 // 397 t = new double[n-1+1, n-1+1]; 398 if( isupperb ) 399 { 400 for(i=0; i<=n-1; i++) 401 { 402 for(i_=i; i_<=n-1;i_++) 403 { 404 t[i,i_] = b[i,i_]; 405 } 406 } 407 } 408 else 409 { 410 for(i=0; i<=n-1; i++) 411 { 412 for(i_=i; i_<=n-1;i_++) 413 { 414 t[i,i_] = b[i_,i]; 415 } 416 } 417 } 418 if( !trfac.spdmatrixcholesky(ref t, n, true) ) 419 { 420 result = false; 421 return result; 422 } 423 424 // 425 // Build U * A * U' in R 426 // 427 w1 = new double[n+1]; 428 w2 = new double[n+1]; 429 w3 = new double[n+1]; 430 r = new double[n-1+1, n-1+1]; 431 for(j=1; j<=n; j++) 432 { 433 434 // 435 // Form w2 = A * u'(j) (here u'(j) is j-th column of U') 436 // 437 i1_ = (j-1) - (1); 438 for(i_=1; i_<=n-j+1;i_++) 439 { 440 w1[i_] = t[j-1,i_+i1_]; 441 } 442 sblas.symmetricmatrixvectormultiply(ref a, isuppera, j-1, n-1, ref w1, 1.0, ref w3); 443 i1_ = (1) - (j); 444 for(i_=j; i_<=n;i_++) 445 { 446 w2[i_] = w3[i_+i1_]; 447 } 448 i1_ = (j-1) - (j); 449 for(i_=j; i_<=n;i_++) 450 { 451 w1[i_] = t[j-1,i_+i1_]; 452 } 453 if( isuppera ) 454 { 455 blas.matrixvectormultiply(ref a, 0, j-2, j-1, n-1, false, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0); 456 } 457 else 458 { 459 blas.matrixvectormultiply(ref a, j-1, n-1, 0, j-2, true, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0); 460 } 461 462 // 463 // Form u(i)*w2 (here u(i) is i-th row of U) 464 // 465 for(i=1; i<=n; i++) 466 { 467 i1_ = (i)-(i-1); 468 v = 0.0; 469 for(i_=i-1; i_<=n-1;i_++) 470 { 471 v += t[i-1,i_]*w2[i_+i1_]; 472 } 473 r[i-1,j-1] = v; 474 } 475 } 476 477 // 478 // Copy R to A 479 // 480 for(i=0; i<=n-1; i++) 481 { 482 for(i_=0; i_<=n-1;i_++) 483 { 484 a[i,i_] = r[i,i_]; 485 } 486 } 487 if( problemtype==2 ) 488 { 489 490 // 491 // Invert U in T 492 // 493 matinv.rmatrixtrinverse(ref t, n, true, false, ref info, ref rep); 494 if( info<=0 ) 495 { 496 result = false; 497 return result; 498 } 499 500 // 501 // Copy U^-1 from T to R 502 // 503 isupperr = true; 504 for(i=0; i<=n-1; i++) 505 { 506 for(j=0; j<=i-1; j++) 507 { 508 r[i,j] = 0; 509 } 510 } 511 for(i=0; i<=n-1; i++) 512 { 513 for(i_=i; i_<=n-1;i_++) 514 { 515 r[i,i_] = t[i,i_]; 516 } 517 } 518 } 519 else 520 { 521 522 // 523 // Copy U from T to R and transpose 524 // 525 isupperr = false; 526 for(i=0; i<=n-1; i++) 527 { 528 for(j=i+1; j<=n-1; j++) 529 { 530 r[i,j] = 0; 531 } 532 } 533 for(i=0; i<=n-1; i++) 534 { 535 for(i_=i; i_<=n-1;i_++) 536 { 537 r[i_,i] = t[i,i_]; 538 } 539 } 540 } 541 } 542 return result; 543 } 436 } 437 438 // 439 // Copy U^-1 from T to R 440 // 441 isupperr = true; 442 for (i = 0; i <= n - 1; i++) { 443 for (j = 0; j <= i - 1; j++) { 444 r[i, j] = 0; 445 } 446 } 447 for (i = 0; i <= n - 1; i++) { 448 for (i_ = i; i_ <= n - 1; i_++) { 449 r[i, i_] = t[i, i_]; 450 } 451 } 452 } else { 453 454 // 455 // Copy U from T to R and transpose 456 // 457 isupperr = false; 458 for (i = 0; i <= n - 1; i++) { 459 for (j = i + 1; j <= n - 1; j++) { 460 r[i, j] = 0; 461 } 462 } 463 for (i = 0; i <= n - 1; i++) { 464 for (i_ = i; i_ <= n - 1; i_++) { 465 r[i_, i] = t[i, i_]; 466 } 467 } 468 } 469 } 470 return result; 544 471 } 472 } 545 473 }
Note: See TracChangeset
for help on using the changeset viewer.
