[2563] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
|
---|
| 3 |
|
---|
| 4 | >>> SOURCE LICENSE >>>
|
---|
| 5 | This program is free software; you can redistribute it and/or modify
|
---|
| 6 | it under the terms of the GNU General Public License as published by
|
---|
| 7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 8 | License, or (at your option) any later version.
|
---|
| 9 |
|
---|
| 10 | This program is distributed in the hope that it will be useful,
|
---|
| 11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 13 | GNU General Public License for more details.
|
---|
| 14 |
|
---|
| 15 | A copy of the GNU General Public License is available at
|
---|
| 16 | http://www.fsf.org/licensing/licenses
|
---|
| 17 |
|
---|
| 18 | >>> END OF LICENSE >>>
|
---|
| 19 | *************************************************************************/
|
---|
| 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 = sevd.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 | }
|
---|
| 174 | 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 | int i_ = 0;
|
---|
| 248 | int i1_ = 0;
|
---|
| 249 |
|
---|
| 250 | System.Diagnostics.Debug.Assert(n>0, "SMatrixGEVDReduce: N<=0!");
|
---|
| 251 | System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "SMatrixGEVDReduce: incorrect ProblemType!");
|
---|
| 252 | result = true;
|
---|
| 253 |
|
---|
| 254 | //
|
---|
| 255 | // Problem 1: A*x = lambda*B*x
|
---|
| 256 | //
|
---|
| 257 | // Reducing to:
|
---|
| 258 | // C*y = lambda*y
|
---|
| 259 | // C = L^(-1) * A * L^(-T)
|
---|
| 260 | // x = L^(-T) * y
|
---|
| 261 | //
|
---|
| 262 | if( problemtype==1 )
|
---|
| 263 | {
|
---|
| 264 |
|
---|
| 265 | //
|
---|
| 266 | // Factorize B in T: B = LL'
|
---|
| 267 | //
|
---|
| 268 | t = new double[n-1+1, n-1+1];
|
---|
| 269 | if( isupperb )
|
---|
| 270 | {
|
---|
| 271 | for(i=0; i<=n-1; i++)
|
---|
| 272 | {
|
---|
| 273 | for(i_=i; i_<=n-1;i_++)
|
---|
| 274 | {
|
---|
| 275 | t[i_,i] = b[i,i_];
|
---|
| 276 | }
|
---|
| 277 | }
|
---|
| 278 | }
|
---|
| 279 | else
|
---|
| 280 | {
|
---|
| 281 | for(i=0; i<=n-1; i++)
|
---|
| 282 | {
|
---|
| 283 | for(i_=0; i_<=i;i_++)
|
---|
| 284 | {
|
---|
| 285 | t[i,i_] = b[i,i_];
|
---|
| 286 | }
|
---|
| 287 | }
|
---|
| 288 | }
|
---|
| 289 | if( !cholesky.spdmatrixcholesky(ref t, n, false) )
|
---|
| 290 | {
|
---|
| 291 | result = false;
|
---|
| 292 | return result;
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | //
|
---|
| 296 | // Invert L in T
|
---|
| 297 | //
|
---|
| 298 | if( !trinverse.rmatrixtrinverse(ref t, n, false, false) )
|
---|
| 299 | {
|
---|
| 300 | result = false;
|
---|
| 301 | return result;
|
---|
| 302 | }
|
---|
| 303 |
|
---|
| 304 | //
|
---|
| 305 | // Build L^(-1) * A * L^(-T) in R
|
---|
| 306 | //
|
---|
| 307 | w1 = new double[n+1];
|
---|
| 308 | w2 = new double[n+1];
|
---|
| 309 | r = new double[n-1+1, n-1+1];
|
---|
| 310 | for(j=1; j<=n; j++)
|
---|
| 311 | {
|
---|
| 312 |
|
---|
| 313 | //
|
---|
| 314 | // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
|
---|
| 315 | //
|
---|
| 316 | i1_ = (0) - (1);
|
---|
| 317 | for(i_=1; i_<=j;i_++)
|
---|
| 318 | {
|
---|
| 319 | w1[i_] = t[j-1,i_+i1_];
|
---|
| 320 | }
|
---|
| 321 | sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j-1, ref w1, 1.0, ref w2);
|
---|
| 322 | if( isuppera )
|
---|
| 323 | {
|
---|
| 324 | blas.matrixvectormultiply(ref a, 0, j-1, j, n-1, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
|
---|
| 325 | }
|
---|
| 326 | else
|
---|
| 327 | {
|
---|
| 328 | blas.matrixvectormultiply(ref a, j, n-1, 0, j-1, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
|
---|
| 329 | }
|
---|
| 330 |
|
---|
| 331 | //
|
---|
| 332 | // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
|
---|
| 333 | //
|
---|
| 334 | for(i=1; i<=n; i++)
|
---|
| 335 | {
|
---|
| 336 | i1_ = (1)-(0);
|
---|
| 337 | v = 0.0;
|
---|
| 338 | for(i_=0; i_<=i-1;i_++)
|
---|
| 339 | {
|
---|
| 340 | v += t[i-1,i_]*w2[i_+i1_];
|
---|
| 341 | }
|
---|
| 342 | r[i-1,j-1] = v;
|
---|
| 343 | }
|
---|
| 344 | }
|
---|
| 345 |
|
---|
| 346 | //
|
---|
| 347 | // Copy R to A
|
---|
| 348 | //
|
---|
| 349 | for(i=0; i<=n-1; i++)
|
---|
| 350 | {
|
---|
| 351 | for(i_=0; i_<=n-1;i_++)
|
---|
| 352 | {
|
---|
| 353 | a[i,i_] = r[i,i_];
|
---|
| 354 | }
|
---|
| 355 | }
|
---|
| 356 |
|
---|
| 357 | //
|
---|
| 358 | // Copy L^(-1) from T to R and transpose
|
---|
| 359 | //
|
---|
| 360 | isupperr = true;
|
---|
| 361 | for(i=0; i<=n-1; i++)
|
---|
| 362 | {
|
---|
| 363 | for(j=0; j<=i-1; j++)
|
---|
| 364 | {
|
---|
| 365 | r[i,j] = 0;
|
---|
| 366 | }
|
---|
| 367 | }
|
---|
| 368 | for(i=0; i<=n-1; i++)
|
---|
| 369 | {
|
---|
| 370 | for(i_=i; i_<=n-1;i_++)
|
---|
| 371 | {
|
---|
| 372 | r[i,i_] = t[i_,i];
|
---|
| 373 | }
|
---|
| 374 | }
|
---|
| 375 | return result;
|
---|
| 376 | }
|
---|
| 377 |
|
---|
| 378 | //
|
---|
| 379 | // Problem 2: A*B*x = lambda*x
|
---|
| 380 | // or
|
---|
| 381 | // problem 3: B*A*x = lambda*x
|
---|
| 382 | //
|
---|
| 383 | // Reducing to:
|
---|
| 384 | // C*y = lambda*y
|
---|
| 385 | // C = U * A * U'
|
---|
| 386 | // B = U'* U
|
---|
| 387 | //
|
---|
| 388 | if( problemtype==2 | problemtype==3 )
|
---|
| 389 | {
|
---|
| 390 |
|
---|
| 391 | //
|
---|
| 392 | // Factorize B in T: B = U'*U
|
---|
| 393 | //
|
---|
| 394 | t = new double[n-1+1, n-1+1];
|
---|
| 395 | if( isupperb )
|
---|
| 396 | {
|
---|
| 397 | for(i=0; i<=n-1; i++)
|
---|
| 398 | {
|
---|
| 399 | for(i_=i; i_<=n-1;i_++)
|
---|
| 400 | {
|
---|
| 401 | t[i,i_] = b[i,i_];
|
---|
| 402 | }
|
---|
| 403 | }
|
---|
| 404 | }
|
---|
| 405 | else
|
---|
| 406 | {
|
---|
| 407 | for(i=0; i<=n-1; i++)
|
---|
| 408 | {
|
---|
| 409 | for(i_=i; i_<=n-1;i_++)
|
---|
| 410 | {
|
---|
| 411 | t[i,i_] = b[i_,i];
|
---|
| 412 | }
|
---|
| 413 | }
|
---|
| 414 | }
|
---|
| 415 | if( !cholesky.spdmatrixcholesky(ref t, n, true) )
|
---|
| 416 | {
|
---|
| 417 | result = false;
|
---|
| 418 | return result;
|
---|
| 419 | }
|
---|
| 420 |
|
---|
| 421 | //
|
---|
| 422 | // Build U * A * U' in R
|
---|
| 423 | //
|
---|
| 424 | w1 = new double[n+1];
|
---|
| 425 | w2 = new double[n+1];
|
---|
| 426 | w3 = new double[n+1];
|
---|
| 427 | r = new double[n-1+1, n-1+1];
|
---|
| 428 | for(j=1; j<=n; j++)
|
---|
| 429 | {
|
---|
| 430 |
|
---|
| 431 | //
|
---|
| 432 | // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
|
---|
| 433 | //
|
---|
| 434 | i1_ = (j-1) - (1);
|
---|
| 435 | for(i_=1; i_<=n-j+1;i_++)
|
---|
| 436 | {
|
---|
| 437 | w1[i_] = t[j-1,i_+i1_];
|
---|
| 438 | }
|
---|
| 439 | sblas.symmetricmatrixvectormultiply(ref a, isuppera, j-1, n-1, ref w1, 1.0, ref w3);
|
---|
| 440 | i1_ = (1) - (j);
|
---|
| 441 | for(i_=j; i_<=n;i_++)
|
---|
| 442 | {
|
---|
| 443 | w2[i_] = w3[i_+i1_];
|
---|
| 444 | }
|
---|
| 445 | i1_ = (j-1) - (j);
|
---|
| 446 | for(i_=j; i_<=n;i_++)
|
---|
| 447 | {
|
---|
| 448 | w1[i_] = t[j-1,i_+i1_];
|
---|
| 449 | }
|
---|
| 450 | if( isuppera )
|
---|
| 451 | {
|
---|
| 452 | 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);
|
---|
| 453 | }
|
---|
| 454 | else
|
---|
| 455 | {
|
---|
| 456 | 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);
|
---|
| 457 | }
|
---|
| 458 |
|
---|
| 459 | //
|
---|
| 460 | // Form u(i)*w2 (here u(i) is i-th row of U)
|
---|
| 461 | //
|
---|
| 462 | for(i=1; i<=n; i++)
|
---|
| 463 | {
|
---|
| 464 | i1_ = (i)-(i-1);
|
---|
| 465 | v = 0.0;
|
---|
| 466 | for(i_=i-1; i_<=n-1;i_++)
|
---|
| 467 | {
|
---|
| 468 | v += t[i-1,i_]*w2[i_+i1_];
|
---|
| 469 | }
|
---|
| 470 | r[i-1,j-1] = v;
|
---|
| 471 | }
|
---|
| 472 | }
|
---|
| 473 |
|
---|
| 474 | //
|
---|
| 475 | // Copy R to A
|
---|
| 476 | //
|
---|
| 477 | for(i=0; i<=n-1; i++)
|
---|
| 478 | {
|
---|
| 479 | for(i_=0; i_<=n-1;i_++)
|
---|
| 480 | {
|
---|
| 481 | a[i,i_] = r[i,i_];
|
---|
| 482 | }
|
---|
| 483 | }
|
---|
| 484 | if( problemtype==2 )
|
---|
| 485 | {
|
---|
| 486 |
|
---|
| 487 | //
|
---|
| 488 | // Invert U in T
|
---|
| 489 | //
|
---|
| 490 | if( !trinverse.rmatrixtrinverse(ref t, n, true, false) )
|
---|
| 491 | {
|
---|
| 492 | result = false;
|
---|
| 493 | return result;
|
---|
| 494 | }
|
---|
| 495 |
|
---|
| 496 | //
|
---|
| 497 | // Copy U^-1 from T to R
|
---|
| 498 | //
|
---|
| 499 | isupperr = true;
|
---|
| 500 | for(i=0; i<=n-1; i++)
|
---|
| 501 | {
|
---|
| 502 | for(j=0; j<=i-1; j++)
|
---|
| 503 | {
|
---|
| 504 | r[i,j] = 0;
|
---|
| 505 | }
|
---|
| 506 | }
|
---|
| 507 | for(i=0; i<=n-1; i++)
|
---|
| 508 | {
|
---|
| 509 | for(i_=i; i_<=n-1;i_++)
|
---|
| 510 | {
|
---|
| 511 | r[i,i_] = t[i,i_];
|
---|
| 512 | }
|
---|
| 513 | }
|
---|
| 514 | }
|
---|
| 515 | else
|
---|
| 516 | {
|
---|
| 517 |
|
---|
| 518 | //
|
---|
| 519 | // Copy U from T to R and transpose
|
---|
| 520 | //
|
---|
| 521 | isupperr = false;
|
---|
| 522 | for(i=0; i<=n-1; i++)
|
---|
| 523 | {
|
---|
| 524 | for(j=i+1; j<=n-1; j++)
|
---|
| 525 | {
|
---|
| 526 | r[i,j] = 0;
|
---|
| 527 | }
|
---|
| 528 | }
|
---|
| 529 | for(i=0; i<=n-1; i++)
|
---|
| 530 | {
|
---|
| 531 | for(i_=i; i_<=n-1;i_++)
|
---|
| 532 | {
|
---|
| 533 | r[i_,i] = t[i,i_];
|
---|
| 534 | }
|
---|
| 535 | }
|
---|
| 536 | }
|
---|
| 537 | }
|
---|
| 538 | return result;
|
---|
| 539 | }
|
---|
| 540 |
|
---|
| 541 |
|
---|
| 542 | public static bool generalizedsymmetricdefiniteevd(double[,] a,
|
---|
| 543 | int n,
|
---|
| 544 | bool isuppera,
|
---|
| 545 | ref double[,] b,
|
---|
| 546 | bool isupperb,
|
---|
| 547 | int zneeded,
|
---|
| 548 | int problemtype,
|
---|
| 549 | ref double[] d,
|
---|
| 550 | ref double[,] z)
|
---|
| 551 | {
|
---|
| 552 | bool result = new bool();
|
---|
| 553 | double[,] r = new double[0,0];
|
---|
| 554 | double[,] t = new double[0,0];
|
---|
| 555 | bool isupperr = new bool();
|
---|
| 556 | int j1 = 0;
|
---|
| 557 | int j2 = 0;
|
---|
| 558 | int j1inc = 0;
|
---|
| 559 | int j2inc = 0;
|
---|
| 560 | int i = 0;
|
---|
| 561 | int j = 0;
|
---|
| 562 | double v = 0;
|
---|
| 563 | int i_ = 0;
|
---|
| 564 |
|
---|
| 565 | a = (double[,])a.Clone();
|
---|
| 566 |
|
---|
| 567 |
|
---|
| 568 | //
|
---|
| 569 | // Reduce and solve
|
---|
| 570 | //
|
---|
| 571 | result = generalizedsymmetricdefiniteevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr);
|
---|
| 572 | if( !result )
|
---|
| 573 | {
|
---|
| 574 | return result;
|
---|
| 575 | }
|
---|
| 576 | result = sevd.symmetricevd(a, n, zneeded, isuppera, ref d, ref t);
|
---|
| 577 | if( !result )
|
---|
| 578 | {
|
---|
| 579 | return result;
|
---|
| 580 | }
|
---|
| 581 |
|
---|
| 582 | //
|
---|
| 583 | // Transform eigenvectors if needed
|
---|
| 584 | //
|
---|
| 585 | if( zneeded!=0 )
|
---|
| 586 | {
|
---|
| 587 |
|
---|
| 588 | //
|
---|
| 589 | // fill Z with zeros
|
---|
| 590 | //
|
---|
| 591 | z = new double[n+1, n+1];
|
---|
| 592 | for(j=1; j<=n; j++)
|
---|
| 593 | {
|
---|
| 594 | z[1,j] = 0.0;
|
---|
| 595 | }
|
---|
| 596 | for(i=2; i<=n; i++)
|
---|
| 597 | {
|
---|
| 598 | for(i_=1; i_<=n;i_++)
|
---|
| 599 | {
|
---|
| 600 | z[i,i_] = z[1,i_];
|
---|
| 601 | }
|
---|
| 602 | }
|
---|
| 603 |
|
---|
| 604 | //
|
---|
| 605 | // Setup R properties
|
---|
| 606 | //
|
---|
| 607 | if( isupperr )
|
---|
| 608 | {
|
---|
| 609 | j1 = 1;
|
---|
| 610 | j2 = n;
|
---|
| 611 | j1inc = +1;
|
---|
| 612 | j2inc = 0;
|
---|
| 613 | }
|
---|
| 614 | else
|
---|
| 615 | {
|
---|
| 616 | j1 = 1;
|
---|
| 617 | j2 = 1;
|
---|
| 618 | j1inc = 0;
|
---|
| 619 | j2inc = +1;
|
---|
| 620 | }
|
---|
| 621 |
|
---|
| 622 | //
|
---|
| 623 | // Calculate R*Z
|
---|
| 624 | //
|
---|
| 625 | for(i=1; i<=n; i++)
|
---|
| 626 | {
|
---|
| 627 | for(j=j1; j<=j2; j++)
|
---|
| 628 | {
|
---|
| 629 | v = r[i,j];
|
---|
| 630 | for(i_=1; i_<=n;i_++)
|
---|
| 631 | {
|
---|
| 632 | z[i,i_] = z[i,i_] + v*t[j,i_];
|
---|
| 633 | }
|
---|
| 634 | }
|
---|
| 635 | j1 = j1+j1inc;
|
---|
| 636 | j2 = j2+j2inc;
|
---|
| 637 | }
|
---|
| 638 | }
|
---|
| 639 | return result;
|
---|
| 640 | }
|
---|
| 641 |
|
---|
| 642 |
|
---|
| 643 | public static bool generalizedsymmetricdefiniteevdreduce(ref double[,] a,
|
---|
| 644 | int n,
|
---|
| 645 | bool isuppera,
|
---|
| 646 | ref double[,] b,
|
---|
| 647 | bool isupperb,
|
---|
| 648 | int problemtype,
|
---|
| 649 | ref double[,] r,
|
---|
| 650 | ref bool isupperr)
|
---|
| 651 | {
|
---|
| 652 | bool result = new bool();
|
---|
| 653 | double[,] t = new double[0,0];
|
---|
| 654 | double[] w1 = new double[0];
|
---|
| 655 | double[] w2 = new double[0];
|
---|
| 656 | double[] w3 = new double[0];
|
---|
| 657 | int i = 0;
|
---|
| 658 | int j = 0;
|
---|
| 659 | double v = 0;
|
---|
| 660 | int i_ = 0;
|
---|
| 661 | int i1_ = 0;
|
---|
| 662 |
|
---|
| 663 | System.Diagnostics.Debug.Assert(n>0, "GeneralizedSymmetricDefiniteEVDReduce: N<=0!");
|
---|
| 664 | System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "GeneralizedSymmetricDefiniteEVDReduce: incorrect ProblemType!");
|
---|
| 665 | result = true;
|
---|
| 666 |
|
---|
| 667 | //
|
---|
| 668 | // Problem 1: A*x = lambda*B*x
|
---|
| 669 | //
|
---|
| 670 | // Reducing to:
|
---|
| 671 | // C*y = lambda*y
|
---|
| 672 | // C = L^(-1) * A * L^(-T)
|
---|
| 673 | // x = L^(-T) * y
|
---|
| 674 | //
|
---|
| 675 | if( problemtype==1 )
|
---|
| 676 | {
|
---|
| 677 |
|
---|
| 678 | //
|
---|
| 679 | // Factorize B in T: B = LL'
|
---|
| 680 | //
|
---|
| 681 | t = new double[n+1, n+1];
|
---|
| 682 | if( isupperb )
|
---|
| 683 | {
|
---|
| 684 | for(i=1; i<=n; i++)
|
---|
| 685 | {
|
---|
| 686 | for(i_=i; i_<=n;i_++)
|
---|
| 687 | {
|
---|
| 688 | t[i_,i] = b[i,i_];
|
---|
| 689 | }
|
---|
| 690 | }
|
---|
| 691 | }
|
---|
| 692 | else
|
---|
| 693 | {
|
---|
| 694 | for(i=1; i<=n; i++)
|
---|
| 695 | {
|
---|
| 696 | for(i_=1; i_<=i;i_++)
|
---|
| 697 | {
|
---|
| 698 | t[i,i_] = b[i,i_];
|
---|
| 699 | }
|
---|
| 700 | }
|
---|
| 701 | }
|
---|
| 702 | if( !cholesky.choleskydecomposition(ref t, n, false) )
|
---|
| 703 | {
|
---|
| 704 | result = false;
|
---|
| 705 | return result;
|
---|
| 706 | }
|
---|
| 707 |
|
---|
| 708 | //
|
---|
| 709 | // Invert L in T
|
---|
| 710 | //
|
---|
| 711 | if( !trinverse.invtriangular(ref t, n, false, false) )
|
---|
| 712 | {
|
---|
| 713 | result = false;
|
---|
| 714 | return result;
|
---|
| 715 | }
|
---|
| 716 |
|
---|
| 717 | //
|
---|
| 718 | // Build L^(-1) * A * L^(-T) in R
|
---|
| 719 | //
|
---|
| 720 | w1 = new double[n+1];
|
---|
| 721 | w2 = new double[n+1];
|
---|
| 722 | r = new double[n+1, n+1];
|
---|
| 723 | for(j=1; j<=n; j++)
|
---|
| 724 | {
|
---|
| 725 |
|
---|
| 726 | //
|
---|
| 727 | // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
|
---|
| 728 | //
|
---|
| 729 | for(i_=1; i_<=j;i_++)
|
---|
| 730 | {
|
---|
| 731 | w1[i_] = t[j,i_];
|
---|
| 732 | }
|
---|
| 733 | sblas.symmetricmatrixvectormultiply(ref a, isuppera, 1, j, ref w1, 1.0, ref w2);
|
---|
| 734 | if( isuppera )
|
---|
| 735 | {
|
---|
| 736 | blas.matrixvectormultiply(ref a, 1, j, j+1, n, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
|
---|
| 737 | }
|
---|
| 738 | else
|
---|
| 739 | {
|
---|
| 740 | blas.matrixvectormultiply(ref a, j+1, n, 1, j, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
|
---|
| 741 | }
|
---|
| 742 |
|
---|
| 743 | //
|
---|
| 744 | // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
|
---|
| 745 | //
|
---|
| 746 | for(i=1; i<=n; i++)
|
---|
| 747 | {
|
---|
| 748 | v = 0.0;
|
---|
| 749 | for(i_=1; i_<=i;i_++)
|
---|
| 750 | {
|
---|
| 751 | v += t[i,i_]*w2[i_];
|
---|
| 752 | }
|
---|
| 753 | r[i,j] = v;
|
---|
| 754 | }
|
---|
| 755 | }
|
---|
| 756 |
|
---|
| 757 | //
|
---|
| 758 | // Copy R to A
|
---|
| 759 | //
|
---|
| 760 | for(i=1; i<=n; i++)
|
---|
| 761 | {
|
---|
| 762 | for(i_=1; i_<=n;i_++)
|
---|
| 763 | {
|
---|
| 764 | a[i,i_] = r[i,i_];
|
---|
| 765 | }
|
---|
| 766 | }
|
---|
| 767 |
|
---|
| 768 | //
|
---|
| 769 | // Copy L^(-1) from T to R and transpose
|
---|
| 770 | //
|
---|
| 771 | isupperr = true;
|
---|
| 772 | for(i=1; i<=n; i++)
|
---|
| 773 | {
|
---|
| 774 | for(j=1; j<=i-1; j++)
|
---|
| 775 | {
|
---|
| 776 | r[i,j] = 0;
|
---|
| 777 | }
|
---|
| 778 | }
|
---|
| 779 | for(i=1; i<=n; i++)
|
---|
| 780 | {
|
---|
| 781 | for(i_=i; i_<=n;i_++)
|
---|
| 782 | {
|
---|
| 783 | r[i,i_] = t[i_,i];
|
---|
| 784 | }
|
---|
| 785 | }
|
---|
| 786 | return result;
|
---|
| 787 | }
|
---|
| 788 |
|
---|
| 789 | //
|
---|
| 790 | // Problem 2: A*B*x = lambda*x
|
---|
| 791 | // or
|
---|
| 792 | // problem 3: B*A*x = lambda*x
|
---|
| 793 | //
|
---|
| 794 | // Reducing to:
|
---|
| 795 | // C*y = lambda*y
|
---|
| 796 | // C = U * A * U'
|
---|
| 797 | // B = U'* U
|
---|
| 798 | //
|
---|
| 799 | if( problemtype==2 | problemtype==3 )
|
---|
| 800 | {
|
---|
| 801 |
|
---|
| 802 | //
|
---|
| 803 | // Factorize B in T: B = U'*U
|
---|
| 804 | //
|
---|
| 805 | t = new double[n+1, n+1];
|
---|
| 806 | if( isupperb )
|
---|
| 807 | {
|
---|
| 808 | for(i=1; i<=n; i++)
|
---|
| 809 | {
|
---|
| 810 | for(i_=i; i_<=n;i_++)
|
---|
| 811 | {
|
---|
| 812 | t[i,i_] = b[i,i_];
|
---|
| 813 | }
|
---|
| 814 | }
|
---|
| 815 | }
|
---|
| 816 | else
|
---|
| 817 | {
|
---|
| 818 | for(i=1; i<=n; i++)
|
---|
| 819 | {
|
---|
| 820 | for(i_=i; i_<=n;i_++)
|
---|
| 821 | {
|
---|
| 822 | t[i,i_] = b[i_,i];
|
---|
| 823 | }
|
---|
| 824 | }
|
---|
| 825 | }
|
---|
| 826 | if( !cholesky.choleskydecomposition(ref t, n, true) )
|
---|
| 827 | {
|
---|
| 828 | result = false;
|
---|
| 829 | return result;
|
---|
| 830 | }
|
---|
| 831 |
|
---|
| 832 | //
|
---|
| 833 | // Build U * A * U' in R
|
---|
| 834 | //
|
---|
| 835 | w1 = new double[n+1];
|
---|
| 836 | w2 = new double[n+1];
|
---|
| 837 | w3 = new double[n+1];
|
---|
| 838 | r = new double[n+1, n+1];
|
---|
| 839 | for(j=1; j<=n; j++)
|
---|
| 840 | {
|
---|
| 841 |
|
---|
| 842 | //
|
---|
| 843 | // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
|
---|
| 844 | //
|
---|
| 845 | i1_ = (j) - (1);
|
---|
| 846 | for(i_=1; i_<=n-j+1;i_++)
|
---|
| 847 | {
|
---|
| 848 | w1[i_] = t[j,i_+i1_];
|
---|
| 849 | }
|
---|
| 850 | sblas.symmetricmatrixvectormultiply(ref a, isuppera, j, n, ref w1, 1.0, ref w3);
|
---|
| 851 | i1_ = (1) - (j);
|
---|
| 852 | for(i_=j; i_<=n;i_++)
|
---|
| 853 | {
|
---|
| 854 | w2[i_] = w3[i_+i1_];
|
---|
| 855 | }
|
---|
| 856 | for(i_=j; i_<=n;i_++)
|
---|
| 857 | {
|
---|
| 858 | w1[i_] = t[j,i_];
|
---|
| 859 | }
|
---|
| 860 | if( isuppera )
|
---|
| 861 | {
|
---|
| 862 | blas.matrixvectormultiply(ref a, 1, j-1, j, n, false, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
|
---|
| 863 | }
|
---|
| 864 | else
|
---|
| 865 | {
|
---|
| 866 | blas.matrixvectormultiply(ref a, j, n, 1, j-1, true, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
|
---|
| 867 | }
|
---|
| 868 |
|
---|
| 869 | //
|
---|
| 870 | // Form u(i)*w2 (here u(i) is i-th row of U)
|
---|
| 871 | //
|
---|
| 872 | for(i=1; i<=n; i++)
|
---|
| 873 | {
|
---|
| 874 | v = 0.0;
|
---|
| 875 | for(i_=i; i_<=n;i_++)
|
---|
| 876 | {
|
---|
| 877 | v += t[i,i_]*w2[i_];
|
---|
| 878 | }
|
---|
| 879 | r[i,j] = v;
|
---|
| 880 | }
|
---|
| 881 | }
|
---|
| 882 |
|
---|
| 883 | //
|
---|
| 884 | // Copy R to A
|
---|
| 885 | //
|
---|
| 886 | for(i=1; i<=n; i++)
|
---|
| 887 | {
|
---|
| 888 | for(i_=1; i_<=n;i_++)
|
---|
| 889 | {
|
---|
| 890 | a[i,i_] = r[i,i_];
|
---|
| 891 | }
|
---|
| 892 | }
|
---|
| 893 | if( problemtype==2 )
|
---|
| 894 | {
|
---|
| 895 |
|
---|
| 896 | //
|
---|
| 897 | // Invert U in T
|
---|
| 898 | //
|
---|
| 899 | if( !trinverse.invtriangular(ref t, n, true, false) )
|
---|
| 900 | {
|
---|
| 901 | result = false;
|
---|
| 902 | return result;
|
---|
| 903 | }
|
---|
| 904 |
|
---|
| 905 | //
|
---|
| 906 | // Copy U^-1 from T to R
|
---|
| 907 | //
|
---|
| 908 | isupperr = true;
|
---|
| 909 | for(i=1; i<=n; i++)
|
---|
| 910 | {
|
---|
| 911 | for(j=1; j<=i-1; j++)
|
---|
| 912 | {
|
---|
| 913 | r[i,j] = 0;
|
---|
| 914 | }
|
---|
| 915 | }
|
---|
| 916 | for(i=1; i<=n; i++)
|
---|
| 917 | {
|
---|
| 918 | for(i_=i; i_<=n;i_++)
|
---|
| 919 | {
|
---|
| 920 | r[i,i_] = t[i,i_];
|
---|
| 921 | }
|
---|
| 922 | }
|
---|
| 923 | }
|
---|
| 924 | else
|
---|
| 925 | {
|
---|
| 926 |
|
---|
| 927 | //
|
---|
| 928 | // Copy U from T to R and transpose
|
---|
| 929 | //
|
---|
| 930 | isupperr = false;
|
---|
| 931 | for(i=1; i<=n; i++)
|
---|
| 932 | {
|
---|
| 933 | for(j=i+1; j<=n; j++)
|
---|
| 934 | {
|
---|
| 935 | r[i,j] = 0;
|
---|
| 936 | }
|
---|
| 937 | }
|
---|
| 938 | for(i=1; i<=n; i++)
|
---|
| 939 | {
|
---|
| 940 | for(i_=i; i_<=n;i_++)
|
---|
| 941 | {
|
---|
| 942 | r[i_,i] = t[i,i_];
|
---|
| 943 | }
|
---|
| 944 | }
|
---|
| 945 | }
|
---|
| 946 | }
|
---|
| 947 | return result;
|
---|
| 948 | }
|
---|
| 949 | }
|
---|
| 950 | }
|
---|