[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
| 3 |
|
---|
| 4 | Contributors:
|
---|
| 5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
| 6 | pseudocode.
|
---|
| 7 |
|
---|
| 8 | See subroutines comments for additional copyrights.
|
---|
| 9 |
|
---|
| 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 >>>
|
---|
| 25 | *************************************************************************/
|
---|
| 26 |
|
---|
| 27 | using System;
|
---|
| 28 |
|
---|
| 29 | namespace alglib
|
---|
| 30 | {
|
---|
| 31 | public class tridiagonal
|
---|
| 32 | {
|
---|
| 33 | /*************************************************************************
|
---|
| 34 | Reduction of a symmetric matrix which is given by its higher or lower
|
---|
| 35 | triangular part to a tridiagonal matrix using orthogonal similarity
|
---|
| 36 | transformation: Q'*A*Q=T.
|
---|
| 37 |
|
---|
| 38 | Input parameters:
|
---|
| 39 | A - matrix to be transformed
|
---|
| 40 | array with elements [0..N-1, 0..N-1].
|
---|
| 41 | N - size of matrix A.
|
---|
| 42 | IsUpper - storage format. If IsUpper = True, then matrix A is given
|
---|
| 43 | by its upper triangle, and the lower triangle is not used
|
---|
| 44 | and not modified by the algorithm, and vice versa
|
---|
| 45 | if IsUpper = False.
|
---|
| 46 |
|
---|
| 47 | Output parameters:
|
---|
| 48 | A - matrices T and Q in compact form (see lower)
|
---|
| 49 | Tau - array of factors which are forming matrices H(i)
|
---|
| 50 | array with elements [0..N-2].
|
---|
| 51 | D - main diagonal of symmetric matrix T.
|
---|
| 52 | array with elements [0..N-1].
|
---|
| 53 | E - secondary diagonal of symmetric matrix T.
|
---|
| 54 | array with elements [0..N-2].
|
---|
| 55 |
|
---|
| 56 |
|
---|
| 57 | If IsUpper=True, the matrix Q is represented as a product of elementary
|
---|
| 58 | reflectors
|
---|
| 59 |
|
---|
| 60 | Q = H(n-2) . . . H(2) H(0).
|
---|
| 61 |
|
---|
| 62 | Each H(i) has the form
|
---|
| 63 |
|
---|
| 64 | H(i) = I - tau * v * v'
|
---|
| 65 |
|
---|
| 66 | where tau is a real scalar, and v is a real vector with
|
---|
| 67 | v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in
|
---|
| 68 | A(0:i-1,i+1), and tau in TAU(i).
|
---|
| 69 |
|
---|
| 70 | If IsUpper=False, the matrix Q is represented as a product of elementary
|
---|
| 71 | reflectors
|
---|
| 72 |
|
---|
| 73 | Q = H(0) H(2) . . . H(n-2).
|
---|
| 74 |
|
---|
| 75 | Each H(i) has the form
|
---|
| 76 |
|
---|
| 77 | H(i) = I - tau * v * v'
|
---|
| 78 |
|
---|
| 79 | where tau is a real scalar, and v is a real vector with
|
---|
| 80 | v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i),
|
---|
| 81 | and tau in TAU(i).
|
---|
| 82 |
|
---|
| 83 | The contents of A on exit are illustrated by the following examples
|
---|
| 84 | with n = 5:
|
---|
| 85 |
|
---|
| 86 | if UPLO = 'U': if UPLO = 'L':
|
---|
| 87 |
|
---|
| 88 | ( d e v1 v2 v3 ) ( d )
|
---|
| 89 | ( d e v2 v3 ) ( e d )
|
---|
| 90 | ( d e v3 ) ( v0 e d )
|
---|
| 91 | ( d e ) ( v0 v1 e d )
|
---|
| 92 | ( d ) ( v0 v1 v2 e d )
|
---|
| 93 |
|
---|
| 94 | where d and e denote diagonal and off-diagonal elements of T, and vi
|
---|
| 95 | denotes an element of the vector defining H(i).
|
---|
| 96 |
|
---|
| 97 | -- LAPACK routine (version 3.0) --
|
---|
| 98 | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
---|
| 99 | Courant Institute, Argonne National Lab, and Rice University
|
---|
| 100 | October 31, 1992
|
---|
| 101 | *************************************************************************/
|
---|
| 102 | public static void smatrixtd(ref double[,] a,
|
---|
| 103 | int n,
|
---|
| 104 | bool isupper,
|
---|
| 105 | ref double[] tau,
|
---|
| 106 | ref double[] d,
|
---|
| 107 | ref double[] e)
|
---|
| 108 | {
|
---|
| 109 | int i = 0;
|
---|
| 110 | double alpha = 0;
|
---|
| 111 | double taui = 0;
|
---|
| 112 | double v = 0;
|
---|
| 113 | double[] t = new double[0];
|
---|
| 114 | double[] t2 = new double[0];
|
---|
| 115 | double[] t3 = new double[0];
|
---|
| 116 | int i_ = 0;
|
---|
| 117 | int i1_ = 0;
|
---|
| 118 |
|
---|
| 119 | if( n<=0 )
|
---|
| 120 | {
|
---|
| 121 | return;
|
---|
| 122 | }
|
---|
| 123 | t = new double[n+1];
|
---|
| 124 | t2 = new double[n+1];
|
---|
| 125 | t3 = new double[n+1];
|
---|
| 126 | if( n>1 )
|
---|
| 127 | {
|
---|
| 128 | tau = new double[n-2+1];
|
---|
| 129 | }
|
---|
| 130 | d = new double[n-1+1];
|
---|
| 131 | if( n>1 )
|
---|
| 132 | {
|
---|
| 133 | e = new double[n-2+1];
|
---|
| 134 | }
|
---|
| 135 | if( isupper )
|
---|
| 136 | {
|
---|
| 137 |
|
---|
| 138 | //
|
---|
| 139 | // Reduce the upper triangle of A
|
---|
| 140 | //
|
---|
| 141 | for(i=n-2; i>=0; i--)
|
---|
| 142 | {
|
---|
| 143 |
|
---|
| 144 | //
|
---|
| 145 | // Generate elementary reflector H() = E - tau * v * v'
|
---|
| 146 | //
|
---|
| 147 | if( i>=1 )
|
---|
| 148 | {
|
---|
| 149 | i1_ = (0) - (2);
|
---|
| 150 | for(i_=2; i_<=i+1;i_++)
|
---|
| 151 | {
|
---|
| 152 | t[i_] = a[i_+i1_,i+1];
|
---|
| 153 | }
|
---|
| 154 | }
|
---|
| 155 | t[1] = a[i,i+1];
|
---|
| 156 | reflections.generatereflection(ref t, i+1, ref taui);
|
---|
| 157 | if( i>=1 )
|
---|
| 158 | {
|
---|
| 159 | i1_ = (2) - (0);
|
---|
| 160 | for(i_=0; i_<=i-1;i_++)
|
---|
| 161 | {
|
---|
| 162 | a[i_,i+1] = t[i_+i1_];
|
---|
| 163 | }
|
---|
| 164 | }
|
---|
| 165 | a[i,i+1] = t[1];
|
---|
| 166 | e[i] = a[i,i+1];
|
---|
| 167 | if( (double)(taui)!=(double)(0) )
|
---|
| 168 | {
|
---|
| 169 |
|
---|
| 170 | //
|
---|
| 171 | // Apply H from both sides to A
|
---|
| 172 | //
|
---|
| 173 | a[i,i+1] = 1;
|
---|
| 174 |
|
---|
| 175 | //
|
---|
| 176 | // Compute x := tau * A * v storing x in TAU
|
---|
| 177 | //
|
---|
| 178 | i1_ = (0) - (1);
|
---|
| 179 | for(i_=1; i_<=i+1;i_++)
|
---|
| 180 | {
|
---|
| 181 | t[i_] = a[i_+i1_,i+1];
|
---|
| 182 | }
|
---|
| 183 | sblas.symmetricmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t3);
|
---|
| 184 | i1_ = (1) - (0);
|
---|
| 185 | for(i_=0; i_<=i;i_++)
|
---|
| 186 | {
|
---|
| 187 | tau[i_] = t3[i_+i1_];
|
---|
| 188 | }
|
---|
| 189 |
|
---|
| 190 | //
|
---|
| 191 | // Compute w := x - 1/2 * tau * (x'*v) * v
|
---|
| 192 | //
|
---|
| 193 | v = 0.0;
|
---|
| 194 | for(i_=0; i_<=i;i_++)
|
---|
| 195 | {
|
---|
| 196 | v += tau[i_]*a[i_,i+1];
|
---|
| 197 | }
|
---|
| 198 | alpha = -(0.5*taui*v);
|
---|
| 199 | for(i_=0; i_<=i;i_++)
|
---|
| 200 | {
|
---|
| 201 | tau[i_] = tau[i_] + alpha*a[i_,i+1];
|
---|
| 202 | }
|
---|
| 203 |
|
---|
| 204 | //
|
---|
| 205 | // Apply the transformation as a rank-2 update:
|
---|
| 206 | // A := A - v * w' - w * v'
|
---|
| 207 | //
|
---|
| 208 | i1_ = (0) - (1);
|
---|
| 209 | for(i_=1; i_<=i+1;i_++)
|
---|
| 210 | {
|
---|
| 211 | t[i_] = a[i_+i1_,i+1];
|
---|
| 212 | }
|
---|
| 213 | i1_ = (0) - (1);
|
---|
| 214 | for(i_=1; i_<=i+1;i_++)
|
---|
| 215 | {
|
---|
| 216 | t3[i_] = tau[i_+i1_];
|
---|
| 217 | }
|
---|
| 218 | sblas.symmetricrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1);
|
---|
| 219 | a[i,i+1] = e[i];
|
---|
| 220 | }
|
---|
| 221 | d[i+1] = a[i+1,i+1];
|
---|
| 222 | tau[i] = taui;
|
---|
| 223 | }
|
---|
| 224 | d[0] = a[0,0];
|
---|
| 225 | }
|
---|
| 226 | else
|
---|
| 227 | {
|
---|
| 228 |
|
---|
| 229 | //
|
---|
| 230 | // Reduce the lower triangle of A
|
---|
| 231 | //
|
---|
| 232 | for(i=0; i<=n-2; i++)
|
---|
| 233 | {
|
---|
| 234 |
|
---|
| 235 | //
|
---|
| 236 | // Generate elementary reflector H = E - tau * v * v'
|
---|
| 237 | //
|
---|
| 238 | i1_ = (i+1) - (1);
|
---|
| 239 | for(i_=1; i_<=n-i-1;i_++)
|
---|
| 240 | {
|
---|
| 241 | t[i_] = a[i_+i1_,i];
|
---|
| 242 | }
|
---|
| 243 | reflections.generatereflection(ref t, n-i-1, ref taui);
|
---|
| 244 | i1_ = (1) - (i+1);
|
---|
| 245 | for(i_=i+1; i_<=n-1;i_++)
|
---|
| 246 | {
|
---|
| 247 | a[i_,i] = t[i_+i1_];
|
---|
| 248 | }
|
---|
| 249 | e[i] = a[i+1,i];
|
---|
| 250 | if( (double)(taui)!=(double)(0) )
|
---|
| 251 | {
|
---|
| 252 |
|
---|
| 253 | //
|
---|
| 254 | // Apply H from both sides to A
|
---|
| 255 | //
|
---|
| 256 | a[i+1,i] = 1;
|
---|
| 257 |
|
---|
| 258 | //
|
---|
| 259 | // Compute x := tau * A * v storing y in TAU
|
---|
| 260 | //
|
---|
| 261 | i1_ = (i+1) - (1);
|
---|
| 262 | for(i_=1; i_<=n-i-1;i_++)
|
---|
| 263 | {
|
---|
| 264 | t[i_] = a[i_+i1_,i];
|
---|
| 265 | }
|
---|
| 266 | sblas.symmetricmatrixvectormultiply(ref a, isupper, i+1, n-1, ref t, taui, ref t2);
|
---|
| 267 | i1_ = (1) - (i);
|
---|
| 268 | for(i_=i; i_<=n-2;i_++)
|
---|
| 269 | {
|
---|
| 270 | tau[i_] = t2[i_+i1_];
|
---|
| 271 | }
|
---|
| 272 |
|
---|
| 273 | //
|
---|
| 274 | // Compute w := x - 1/2 * tau * (x'*v) * v
|
---|
| 275 | //
|
---|
| 276 | i1_ = (i+1)-(i);
|
---|
| 277 | v = 0.0;
|
---|
| 278 | for(i_=i; i_<=n-2;i_++)
|
---|
| 279 | {
|
---|
| 280 | v += tau[i_]*a[i_+i1_,i];
|
---|
| 281 | }
|
---|
| 282 | alpha = -(0.5*taui*v);
|
---|
| 283 | i1_ = (i+1) - (i);
|
---|
| 284 | for(i_=i; i_<=n-2;i_++)
|
---|
| 285 | {
|
---|
| 286 | tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
|
---|
| 287 | }
|
---|
| 288 |
|
---|
| 289 | //
|
---|
| 290 | // Apply the transformation as a rank-2 update:
|
---|
| 291 | // A := A - v * w' - w * v'
|
---|
| 292 | //
|
---|
| 293 | //
|
---|
| 294 | i1_ = (i+1) - (1);
|
---|
| 295 | for(i_=1; i_<=n-i-1;i_++)
|
---|
| 296 | {
|
---|
| 297 | t[i_] = a[i_+i1_,i];
|
---|
| 298 | }
|
---|
| 299 | i1_ = (i) - (1);
|
---|
| 300 | for(i_=1; i_<=n-i-1;i_++)
|
---|
| 301 | {
|
---|
| 302 | t2[i_] = tau[i_+i1_];
|
---|
| 303 | }
|
---|
| 304 | sblas.symmetricrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1);
|
---|
| 305 | a[i+1,i] = e[i];
|
---|
| 306 | }
|
---|
| 307 | d[i] = a[i,i];
|
---|
| 308 | tau[i] = taui;
|
---|
| 309 | }
|
---|
| 310 | d[n-1] = a[n-1,n-1];
|
---|
| 311 | }
|
---|
| 312 | }
|
---|
| 313 |
|
---|
| 314 |
|
---|
| 315 | /*************************************************************************
|
---|
| 316 | Unpacking matrix Q which reduces symmetric matrix to a tridiagonal
|
---|
| 317 | form.
|
---|
| 318 |
|
---|
| 319 | Input parameters:
|
---|
| 320 | A - the result of a SMatrixTD subroutine
|
---|
| 321 | N - size of matrix A.
|
---|
| 322 | IsUpper - storage format (a parameter of SMatrixTD subroutine)
|
---|
| 323 | Tau - the result of a SMatrixTD subroutine
|
---|
| 324 |
|
---|
| 325 | Output parameters:
|
---|
| 326 | Q - transformation matrix.
|
---|
| 327 | array with elements [0..N-1, 0..N-1].
|
---|
| 328 |
|
---|
| 329 | -- ALGLIB --
|
---|
| 330 | Copyright 2005-2008 by Bochkanov Sergey
|
---|
| 331 | *************************************************************************/
|
---|
| 332 | public static void smatrixtdunpackq(ref double[,] a,
|
---|
| 333 | int n,
|
---|
| 334 | bool isupper,
|
---|
| 335 | ref double[] tau,
|
---|
| 336 | ref double[,] q)
|
---|
| 337 | {
|
---|
| 338 | int i = 0;
|
---|
| 339 | int j = 0;
|
---|
| 340 | double[] v = new double[0];
|
---|
| 341 | double[] work = new double[0];
|
---|
| 342 | int i_ = 0;
|
---|
| 343 | int i1_ = 0;
|
---|
| 344 |
|
---|
| 345 | if( n==0 )
|
---|
| 346 | {
|
---|
| 347 | return;
|
---|
| 348 | }
|
---|
| 349 |
|
---|
| 350 | //
|
---|
| 351 | // init
|
---|
| 352 | //
|
---|
| 353 | q = new double[n-1+1, n-1+1];
|
---|
| 354 | v = new double[n+1];
|
---|
| 355 | work = new double[n-1+1];
|
---|
| 356 | for(i=0; i<=n-1; i++)
|
---|
| 357 | {
|
---|
| 358 | for(j=0; j<=n-1; j++)
|
---|
| 359 | {
|
---|
| 360 | if( i==j )
|
---|
| 361 | {
|
---|
| 362 | q[i,j] = 1;
|
---|
| 363 | }
|
---|
| 364 | else
|
---|
| 365 | {
|
---|
| 366 | q[i,j] = 0;
|
---|
| 367 | }
|
---|
| 368 | }
|
---|
| 369 | }
|
---|
| 370 |
|
---|
| 371 | //
|
---|
| 372 | // unpack Q
|
---|
| 373 | //
|
---|
| 374 | if( isupper )
|
---|
| 375 | {
|
---|
| 376 | for(i=0; i<=n-2; i++)
|
---|
| 377 | {
|
---|
| 378 |
|
---|
| 379 | //
|
---|
| 380 | // Apply H(i)
|
---|
| 381 | //
|
---|
| 382 | i1_ = (0) - (1);
|
---|
| 383 | for(i_=1; i_<=i+1;i_++)
|
---|
| 384 | {
|
---|
| 385 | v[i_] = a[i_+i1_,i+1];
|
---|
| 386 | }
|
---|
| 387 | v[i+1] = 1;
|
---|
| 388 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, 0, i, 0, n-1, ref work);
|
---|
| 389 | }
|
---|
| 390 | }
|
---|
| 391 | else
|
---|
| 392 | {
|
---|
| 393 | for(i=n-2; i>=0; i--)
|
---|
| 394 | {
|
---|
| 395 |
|
---|
| 396 | //
|
---|
| 397 | // Apply H(i)
|
---|
| 398 | //
|
---|
| 399 | i1_ = (i+1) - (1);
|
---|
| 400 | for(i_=1; i_<=n-i-1;i_++)
|
---|
| 401 | {
|
---|
| 402 | v[i_] = a[i_+i1_,i];
|
---|
| 403 | }
|
---|
| 404 | v[1] = 1;
|
---|
| 405 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work);
|
---|
| 406 | }
|
---|
| 407 | }
|
---|
| 408 | }
|
---|
| 409 |
|
---|
| 410 |
|
---|
| 411 | public static void totridiagonal(ref double[,] a,
|
---|
| 412 | int n,
|
---|
| 413 | bool isupper,
|
---|
| 414 | ref double[] tau,
|
---|
| 415 | ref double[] d,
|
---|
| 416 | ref double[] e)
|
---|
| 417 | {
|
---|
| 418 | int i = 0;
|
---|
| 419 | int ip1 = 0;
|
---|
| 420 | int im1 = 0;
|
---|
| 421 | int nmi = 0;
|
---|
| 422 | int nm1 = 0;
|
---|
| 423 | double alpha = 0;
|
---|
| 424 | double taui = 0;
|
---|
| 425 | double v = 0;
|
---|
| 426 | double[] t = new double[0];
|
---|
| 427 | double[] t2 = new double[0];
|
---|
| 428 | double[] t3 = new double[0];
|
---|
| 429 | int i_ = 0;
|
---|
| 430 | int i1_ = 0;
|
---|
| 431 |
|
---|
| 432 | if( n<=0 )
|
---|
| 433 | {
|
---|
| 434 | return;
|
---|
| 435 | }
|
---|
| 436 | t = new double[n+1];
|
---|
| 437 | t2 = new double[n+1];
|
---|
| 438 | t3 = new double[n+1];
|
---|
| 439 | tau = new double[Math.Max(1, n-1)+1];
|
---|
| 440 | d = new double[n+1];
|
---|
| 441 | e = new double[Math.Max(1, n-1)+1];
|
---|
| 442 | if( isupper )
|
---|
| 443 | {
|
---|
| 444 |
|
---|
| 445 | //
|
---|
| 446 | // Reduce the upper triangle of A
|
---|
| 447 | //
|
---|
| 448 | for(i=n-1; i>=1; i--)
|
---|
| 449 | {
|
---|
| 450 |
|
---|
| 451 | //
|
---|
| 452 | // Generate elementary reflector H(i) = I - tau * v * v'
|
---|
| 453 | // to annihilate A(1:i-1,i+1)
|
---|
| 454 | //
|
---|
| 455 | // DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI );
|
---|
| 456 | //
|
---|
| 457 | ip1 = i+1;
|
---|
| 458 | im1 = i-1;
|
---|
| 459 | if( i>=2 )
|
---|
| 460 | {
|
---|
| 461 | i1_ = (1) - (2);
|
---|
| 462 | for(i_=2; i_<=i;i_++)
|
---|
| 463 | {
|
---|
| 464 | t[i_] = a[i_+i1_,ip1];
|
---|
| 465 | }
|
---|
| 466 | }
|
---|
| 467 | t[1] = a[i,ip1];
|
---|
| 468 | reflections.generatereflection(ref t, i, ref taui);
|
---|
| 469 | if( i>=2 )
|
---|
| 470 | {
|
---|
| 471 | i1_ = (2) - (1);
|
---|
| 472 | for(i_=1; i_<=im1;i_++)
|
---|
| 473 | {
|
---|
| 474 | a[i_,ip1] = t[i_+i1_];
|
---|
| 475 | }
|
---|
| 476 | }
|
---|
| 477 | a[i,ip1] = t[1];
|
---|
| 478 | e[i] = a[i,i+1];
|
---|
| 479 | if( (double)(taui)!=(double)(0) )
|
---|
| 480 | {
|
---|
| 481 |
|
---|
| 482 | //
|
---|
| 483 | // Apply H(i) from both sides to A(1:i,1:i)
|
---|
| 484 | //
|
---|
| 485 | a[i,i+1] = 1;
|
---|
| 486 |
|
---|
| 487 | //
|
---|
| 488 | // Compute x := tau * A * v storing x in TAU(1:i)
|
---|
| 489 | //
|
---|
| 490 | // DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, TAU, 1 );
|
---|
| 491 | //
|
---|
| 492 | ip1 = i+1;
|
---|
| 493 | for(i_=1; i_<=i;i_++)
|
---|
| 494 | {
|
---|
| 495 | t[i_] = a[i_,ip1];
|
---|
| 496 | }
|
---|
| 497 | sblas.symmetricmatrixvectormultiply(ref a, isupper, 1, i, ref t, taui, ref tau);
|
---|
| 498 |
|
---|
| 499 | //
|
---|
| 500 | // Compute w := x - 1/2 * tau * (x'*v) * v
|
---|
| 501 | //
|
---|
| 502 | ip1 = i+1;
|
---|
| 503 | v = 0.0;
|
---|
| 504 | for(i_=1; i_<=i;i_++)
|
---|
| 505 | {
|
---|
| 506 | v += tau[i_]*a[i_,ip1];
|
---|
| 507 | }
|
---|
| 508 | alpha = -(0.5*taui*v);
|
---|
| 509 | for(i_=1; i_<=i;i_++)
|
---|
| 510 | {
|
---|
| 511 | tau[i_] = tau[i_] + alpha*a[i_,ip1];
|
---|
| 512 | }
|
---|
| 513 |
|
---|
| 514 | //
|
---|
| 515 | // Apply the transformation as a rank-2 update:
|
---|
| 516 | // A := A - v * w' - w * v'
|
---|
| 517 | //
|
---|
| 518 | // DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, LDA );
|
---|
| 519 | //
|
---|
| 520 | for(i_=1; i_<=i;i_++)
|
---|
| 521 | {
|
---|
| 522 | t[i_] = a[i_,ip1];
|
---|
| 523 | }
|
---|
| 524 | sblas.symmetricrank2update(ref a, isupper, 1, i, ref t, ref tau, ref t2, -1);
|
---|
| 525 | a[i,i+1] = e[i];
|
---|
| 526 | }
|
---|
| 527 | d[i+1] = a[i+1,i+1];
|
---|
| 528 | tau[i] = taui;
|
---|
| 529 | }
|
---|
| 530 | d[1] = a[1,1];
|
---|
| 531 | }
|
---|
| 532 | else
|
---|
| 533 | {
|
---|
| 534 |
|
---|
| 535 | //
|
---|
| 536 | // Reduce the lower triangle of A
|
---|
| 537 | //
|
---|
| 538 | for(i=1; i<=n-1; i++)
|
---|
| 539 | {
|
---|
| 540 |
|
---|
| 541 | //
|
---|
| 542 | // Generate elementary reflector H(i) = I - tau * v * v'
|
---|
| 543 | // to annihilate A(i+2:n,i)
|
---|
| 544 | //
|
---|
| 545 | //DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, TAUI );
|
---|
| 546 | //
|
---|
| 547 | nmi = n-i;
|
---|
| 548 | ip1 = i+1;
|
---|
| 549 | i1_ = (ip1) - (1);
|
---|
| 550 | for(i_=1; i_<=nmi;i_++)
|
---|
| 551 | {
|
---|
| 552 | t[i_] = a[i_+i1_,i];
|
---|
| 553 | }
|
---|
| 554 | reflections.generatereflection(ref t, nmi, ref taui);
|
---|
| 555 | i1_ = (1) - (ip1);
|
---|
| 556 | for(i_=ip1; i_<=n;i_++)
|
---|
| 557 | {
|
---|
| 558 | a[i_,i] = t[i_+i1_];
|
---|
| 559 | }
|
---|
| 560 | e[i] = a[i+1,i];
|
---|
| 561 | if( (double)(taui)!=(double)(0) )
|
---|
| 562 | {
|
---|
| 563 |
|
---|
| 564 | //
|
---|
| 565 | // Apply H(i) from both sides to A(i+1:n,i+1:n)
|
---|
| 566 | //
|
---|
| 567 | a[i+1,i] = 1;
|
---|
| 568 |
|
---|
| 569 | //
|
---|
| 570 | // Compute x := tau * A * v storing y in TAU(i:n-1)
|
---|
| 571 | //
|
---|
| 572 | //DSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO, TAU( I ), 1 );
|
---|
| 573 | //
|
---|
| 574 | ip1 = i+1;
|
---|
| 575 | nmi = n-i;
|
---|
| 576 | nm1 = n-1;
|
---|
| 577 | i1_ = (ip1) - (1);
|
---|
| 578 | for(i_=1; i_<=nmi;i_++)
|
---|
| 579 | {
|
---|
| 580 | t[i_] = a[i_+i1_,i];
|
---|
| 581 | }
|
---|
| 582 | sblas.symmetricmatrixvectormultiply(ref a, isupper, i+1, n, ref t, taui, ref t2);
|
---|
| 583 | i1_ = (1) - (i);
|
---|
| 584 | for(i_=i; i_<=nm1;i_++)
|
---|
| 585 | {
|
---|
| 586 | tau[i_] = t2[i_+i1_];
|
---|
| 587 | }
|
---|
| 588 |
|
---|
| 589 | //
|
---|
| 590 | // Compute w := x - 1/2 * tau * (x'*v) * v
|
---|
| 591 | //
|
---|
| 592 | nm1 = n-1;
|
---|
| 593 | ip1 = i+1;
|
---|
| 594 | i1_ = (ip1)-(i);
|
---|
| 595 | v = 0.0;
|
---|
| 596 | for(i_=i; i_<=nm1;i_++)
|
---|
| 597 | {
|
---|
| 598 | v += tau[i_]*a[i_+i1_,i];
|
---|
| 599 | }
|
---|
| 600 | alpha = -(0.5*taui*v);
|
---|
| 601 | i1_ = (ip1) - (i);
|
---|
| 602 | for(i_=i; i_<=nm1;i_++)
|
---|
| 603 | {
|
---|
| 604 | tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
|
---|
| 605 | }
|
---|
| 606 |
|
---|
| 607 | //
|
---|
| 608 | // Apply the transformation as a rank-2 update:
|
---|
| 609 | // A := A - v * w' - w * v'
|
---|
| 610 | //
|
---|
| 611 | //DSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, A( I+1, I+1 ), LDA );
|
---|
| 612 | //
|
---|
| 613 | nm1 = n-1;
|
---|
| 614 | nmi = n-i;
|
---|
| 615 | ip1 = i+1;
|
---|
| 616 | i1_ = (ip1) - (1);
|
---|
| 617 | for(i_=1; i_<=nmi;i_++)
|
---|
| 618 | {
|
---|
| 619 | t[i_] = a[i_+i1_,i];
|
---|
| 620 | }
|
---|
| 621 | i1_ = (i) - (1);
|
---|
| 622 | for(i_=1; i_<=nmi;i_++)
|
---|
| 623 | {
|
---|
| 624 | t2[i_] = tau[i_+i1_];
|
---|
| 625 | }
|
---|
| 626 | sblas.symmetricrank2update(ref a, isupper, i+1, n, ref t, ref t2, ref t3, -1);
|
---|
| 627 | a[i+1,i] = e[i];
|
---|
| 628 | }
|
---|
| 629 | d[i] = a[i,i];
|
---|
| 630 | tau[i] = taui;
|
---|
| 631 | }
|
---|
| 632 | d[n] = a[n,n];
|
---|
| 633 | }
|
---|
| 634 | }
|
---|
| 635 |
|
---|
| 636 |
|
---|
| 637 | public static void unpackqfromtridiagonal(ref double[,] a,
|
---|
| 638 | int n,
|
---|
| 639 | bool isupper,
|
---|
| 640 | ref double[] tau,
|
---|
| 641 | ref double[,] q)
|
---|
| 642 | {
|
---|
| 643 | int i = 0;
|
---|
| 644 | int j = 0;
|
---|
| 645 | int ip1 = 0;
|
---|
| 646 | int nmi = 0;
|
---|
| 647 | double[] v = new double[0];
|
---|
| 648 | double[] work = new double[0];
|
---|
| 649 | int i_ = 0;
|
---|
| 650 | int i1_ = 0;
|
---|
| 651 |
|
---|
| 652 | if( n==0 )
|
---|
| 653 | {
|
---|
| 654 | return;
|
---|
| 655 | }
|
---|
| 656 |
|
---|
| 657 | //
|
---|
| 658 | // init
|
---|
| 659 | //
|
---|
| 660 | q = new double[n+1, n+1];
|
---|
| 661 | v = new double[n+1];
|
---|
| 662 | work = new double[n+1];
|
---|
| 663 | for(i=1; i<=n; i++)
|
---|
| 664 | {
|
---|
| 665 | for(j=1; j<=n; j++)
|
---|
| 666 | {
|
---|
| 667 | if( i==j )
|
---|
| 668 | {
|
---|
| 669 | q[i,j] = 1;
|
---|
| 670 | }
|
---|
| 671 | else
|
---|
| 672 | {
|
---|
| 673 | q[i,j] = 0;
|
---|
| 674 | }
|
---|
| 675 | }
|
---|
| 676 | }
|
---|
| 677 |
|
---|
| 678 | //
|
---|
| 679 | // unpack Q
|
---|
| 680 | //
|
---|
| 681 | if( isupper )
|
---|
| 682 | {
|
---|
| 683 | for(i=1; i<=n-1; i++)
|
---|
| 684 | {
|
---|
| 685 |
|
---|
| 686 | //
|
---|
| 687 | // Apply H(i)
|
---|
| 688 | //
|
---|
| 689 | ip1 = i+1;
|
---|
| 690 | for(i_=1; i_<=i;i_++)
|
---|
| 691 | {
|
---|
| 692 | v[i_] = a[i_,ip1];
|
---|
| 693 | }
|
---|
| 694 | v[i] = 1;
|
---|
| 695 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, 1, i, 1, n, ref work);
|
---|
| 696 | }
|
---|
| 697 | }
|
---|
| 698 | else
|
---|
| 699 | {
|
---|
| 700 | for(i=n-1; i>=1; i--)
|
---|
| 701 | {
|
---|
| 702 |
|
---|
| 703 | //
|
---|
| 704 | // Apply H(i)
|
---|
| 705 | //
|
---|
| 706 | ip1 = i+1;
|
---|
| 707 | nmi = n-i;
|
---|
| 708 | i1_ = (ip1) - (1);
|
---|
| 709 | for(i_=1; i_<=nmi;i_++)
|
---|
| 710 | {
|
---|
| 711 | v[i_] = a[i_+i1_,i];
|
---|
| 712 | }
|
---|
| 713 | v[1] = 1;
|
---|
| 714 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n, 1, n, ref work);
|
---|
| 715 | }
|
---|
| 716 | }
|
---|
| 717 | }
|
---|
| 718 | }
|
---|
| 719 | }
|
---|