[2154] | 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 |
|
---|
[2430] | 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.
|
---|
[2154] | 15 |
|
---|
[2430] | 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.
|
---|
[2154] | 20 |
|
---|
[2430] | 21 | A copy of the GNU General Public License is available at
|
---|
| 22 | http://www.fsf.org/licensing/licenses
|
---|
[2154] | 23 |
|
---|
[2430] | 24 | >>> END OF LICENSE >>>
|
---|
[2154] | 25 | *************************************************************************/
|
---|
| 26 |
|
---|
| 27 | using System;
|
---|
| 28 |
|
---|
[2430] | 29 | namespace alglib
|
---|
[2154] | 30 | {
|
---|
[2430] | 31 | public class bdsvd
|
---|
| 32 | {
|
---|
| 33 | /*************************************************************************
|
---|
| 34 | Singular value decomposition of a bidiagonal matrix (extended algorithm)
|
---|
[2154] | 35 |
|
---|
[2430] | 36 | The algorithm performs the singular value decomposition of a bidiagonal
|
---|
| 37 | matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and P -
|
---|
| 38 | orthogonal matrices, S - diagonal matrix with non-negative elements on the
|
---|
| 39 | main diagonal, in descending order.
|
---|
[2154] | 40 |
|
---|
[2430] | 41 | The algorithm finds singular values. In addition, the algorithm can
|
---|
| 42 | calculate matrices Q and P (more precisely, not the matrices, but their
|
---|
| 43 | product with given matrices U and VT - U*Q and (P^T)*VT)). Of course,
|
---|
| 44 | matrices U and VT can be of any type, including identity. Furthermore, the
|
---|
| 45 | algorithm can calculate Q'*C (this product is calculated more effectively
|
---|
| 46 | than U*Q, because this calculation operates with rows instead of matrix
|
---|
| 47 | columns).
|
---|
[2154] | 48 |
|
---|
[2430] | 49 | The feature of the algorithm is its ability to find all singular values
|
---|
| 50 | including those which are arbitrarily close to 0 with relative accuracy
|
---|
| 51 | close to machine precision. If the parameter IsFractionalAccuracyRequired
|
---|
| 52 | is set to True, all singular values will have high relative accuracy close
|
---|
| 53 | to machine precision. If the parameter is set to False, only the biggest
|
---|
| 54 | singular value will have relative accuracy close to machine precision.
|
---|
| 55 | The absolute error of other singular values is equal to the absolute error
|
---|
| 56 | of the biggest singular value.
|
---|
[2154] | 57 |
|
---|
[2430] | 58 | Input parameters:
|
---|
| 59 | D - main diagonal of matrix B.
|
---|
| 60 | Array whose index ranges within [0..N-1].
|
---|
| 61 | E - superdiagonal (or subdiagonal) of matrix B.
|
---|
| 62 | Array whose index ranges within [0..N-2].
|
---|
| 63 | N - size of matrix B.
|
---|
| 64 | IsUpper - True, if the matrix is upper bidiagonal.
|
---|
| 65 | IsFractionalAccuracyRequired -
|
---|
| 66 | accuracy to search singular values with.
|
---|
| 67 | U - matrix to be multiplied by Q.
|
---|
| 68 | Array whose indexes range within [0..NRU-1, 0..N-1].
|
---|
| 69 | The matrix can be bigger, in that case only the submatrix
|
---|
| 70 | [0..NRU-1, 0..N-1] will be multiplied by Q.
|
---|
| 71 | NRU - number of rows in matrix U.
|
---|
| 72 | C - matrix to be multiplied by Q'.
|
---|
| 73 | Array whose indexes range within [0..N-1, 0..NCC-1].
|
---|
| 74 | The matrix can be bigger, in that case only the submatrix
|
---|
| 75 | [0..N-1, 0..NCC-1] will be multiplied by Q'.
|
---|
| 76 | NCC - number of columns in matrix C.
|
---|
| 77 | VT - matrix to be multiplied by P^T.
|
---|
| 78 | Array whose indexes range within [0..N-1, 0..NCVT-1].
|
---|
| 79 | The matrix can be bigger, in that case only the submatrix
|
---|
| 80 | [0..N-1, 0..NCVT-1] will be multiplied by P^T.
|
---|
| 81 | NCVT - number of columns in matrix VT.
|
---|
[2154] | 82 |
|
---|
[2430] | 83 | Output parameters:
|
---|
| 84 | D - singular values of matrix B in descending order.
|
---|
| 85 | U - if NRU>0, contains matrix U*Q.
|
---|
| 86 | VT - if NCVT>0, contains matrix (P^T)*VT.
|
---|
| 87 | C - if NCC>0, contains matrix Q'*C.
|
---|
[2154] | 88 |
|
---|
[2430] | 89 | Result:
|
---|
| 90 | True, if the algorithm has converged.
|
---|
| 91 | False, if the algorithm hasn't converged (rare case).
|
---|
[2154] | 92 |
|
---|
[2430] | 93 | Additional information:
|
---|
| 94 | The type of convergence is controlled by the internal parameter TOL.
|
---|
| 95 | If the parameter is greater than 0, the singular values will have
|
---|
| 96 | relative accuracy TOL. If TOL<0, the singular values will have
|
---|
| 97 | absolute accuracy ABS(TOL)*norm(B).
|
---|
| 98 | By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon,
|
---|
| 99 | where Epsilon is the machine precision. It is not recommended to use
|
---|
| 100 | TOL less than 10*Epsilon since this will considerably slow down the
|
---|
| 101 | algorithm and may not lead to error decreasing.
|
---|
| 102 | History:
|
---|
| 103 | * 31 March, 2007.
|
---|
| 104 | changed MAXITR from 6 to 12.
|
---|
[2154] | 105 |
|
---|
[2430] | 106 | -- LAPACK routine (version 3.0) --
|
---|
| 107 | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
---|
| 108 | Courant Institute, Argonne National Lab, and Rice University
|
---|
| 109 | October 31, 1999.
|
---|
| 110 | *************************************************************************/
|
---|
| 111 | public static bool rmatrixbdsvd(ref double[] d,
|
---|
| 112 | double[] e,
|
---|
| 113 | int n,
|
---|
| 114 | bool isupper,
|
---|
| 115 | bool isfractionalaccuracyrequired,
|
---|
| 116 | ref double[,] u,
|
---|
| 117 | int nru,
|
---|
| 118 | ref double[,] c,
|
---|
| 119 | int ncc,
|
---|
| 120 | ref double[,] vt,
|
---|
| 121 | int ncvt)
|
---|
| 122 | {
|
---|
| 123 | bool result = new bool();
|
---|
| 124 | double[] d1 = new double[0];
|
---|
| 125 | double[] e1 = new double[0];
|
---|
| 126 | int i_ = 0;
|
---|
| 127 | int i1_ = 0;
|
---|
[2154] | 128 |
|
---|
[2430] | 129 | e = (double[])e.Clone();
|
---|
[2154] | 130 |
|
---|
[2430] | 131 | d1 = new double[n+1];
|
---|
[2154] | 132 | i1_ = (0) - (1);
|
---|
[2430] | 133 | for(i_=1; i_<=n;i_++)
|
---|
[2154] | 134 | {
|
---|
[2430] | 135 | d1[i_] = d[i_+i1_];
|
---|
[2154] | 136 | }
|
---|
[2430] | 137 | if( n>1 )
|
---|
| 138 | {
|
---|
| 139 | e1 = new double[n-1+1];
|
---|
| 140 | i1_ = (0) - (1);
|
---|
| 141 | for(i_=1; i_<=n-1;i_++)
|
---|
| 142 | {
|
---|
| 143 | e1[i_] = e[i_+i1_];
|
---|
| 144 | }
|
---|
| 145 | }
|
---|
| 146 | result = bidiagonalsvddecompositioninternal(ref d1, e1, n, isupper, isfractionalaccuracyrequired, ref u, 0, nru, ref c, 0, ncc, ref vt, 0, ncvt);
|
---|
| 147 | i1_ = (1) - (0);
|
---|
| 148 | for(i_=0; i_<=n-1;i_++)
|
---|
| 149 | {
|
---|
| 150 | d[i_] = d1[i_+i1_];
|
---|
| 151 | }
|
---|
| 152 | return result;
|
---|
[2154] | 153 | }
|
---|
| 154 |
|
---|
| 155 |
|
---|
[2430] | 156 | public static bool bidiagonalsvddecomposition(ref double[] d,
|
---|
| 157 | double[] e,
|
---|
| 158 | int n,
|
---|
| 159 | bool isupper,
|
---|
| 160 | bool isfractionalaccuracyrequired,
|
---|
| 161 | ref double[,] u,
|
---|
| 162 | int nru,
|
---|
| 163 | ref double[,] c,
|
---|
| 164 | int ncc,
|
---|
| 165 | ref double[,] vt,
|
---|
| 166 | int ncvt)
|
---|
| 167 | {
|
---|
| 168 | bool result = new bool();
|
---|
[2154] | 169 |
|
---|
[2430] | 170 | e = (double[])e.Clone();
|
---|
[2154] | 171 |
|
---|
[2430] | 172 | result = bidiagonalsvddecompositioninternal(ref d, e, n, isupper, isfractionalaccuracyrequired, ref u, 1, nru, ref c, 1, ncc, ref vt, 1, ncvt);
|
---|
| 173 | return result;
|
---|
| 174 | }
|
---|
[2154] | 175 |
|
---|
| 176 |
|
---|
[2430] | 177 | /*************************************************************************
|
---|
| 178 | Internal working subroutine for bidiagonal decomposition
|
---|
| 179 | *************************************************************************/
|
---|
| 180 | private static bool bidiagonalsvddecompositioninternal(ref double[] d,
|
---|
| 181 | double[] e,
|
---|
| 182 | int n,
|
---|
| 183 | bool isupper,
|
---|
| 184 | bool isfractionalaccuracyrequired,
|
---|
| 185 | ref double[,] u,
|
---|
| 186 | int ustart,
|
---|
| 187 | int nru,
|
---|
| 188 | ref double[,] c,
|
---|
| 189 | int cstart,
|
---|
| 190 | int ncc,
|
---|
| 191 | ref double[,] vt,
|
---|
| 192 | int vstart,
|
---|
| 193 | int ncvt)
|
---|
| 194 | {
|
---|
| 195 | bool result = new bool();
|
---|
| 196 | int i = 0;
|
---|
| 197 | int idir = 0;
|
---|
| 198 | int isub = 0;
|
---|
| 199 | int iter = 0;
|
---|
| 200 | int j = 0;
|
---|
| 201 | int ll = 0;
|
---|
| 202 | int lll = 0;
|
---|
| 203 | int m = 0;
|
---|
| 204 | int maxit = 0;
|
---|
| 205 | int oldll = 0;
|
---|
| 206 | int oldm = 0;
|
---|
| 207 | double abse = 0;
|
---|
| 208 | double abss = 0;
|
---|
| 209 | double cosl = 0;
|
---|
| 210 | double cosr = 0;
|
---|
| 211 | double cs = 0;
|
---|
| 212 | double eps = 0;
|
---|
| 213 | double f = 0;
|
---|
| 214 | double g = 0;
|
---|
| 215 | double h = 0;
|
---|
| 216 | double mu = 0;
|
---|
| 217 | double oldcs = 0;
|
---|
| 218 | double oldsn = 0;
|
---|
| 219 | double r = 0;
|
---|
| 220 | double shift = 0;
|
---|
| 221 | double sigmn = 0;
|
---|
| 222 | double sigmx = 0;
|
---|
| 223 | double sinl = 0;
|
---|
| 224 | double sinr = 0;
|
---|
| 225 | double sll = 0;
|
---|
| 226 | double smax = 0;
|
---|
| 227 | double smin = 0;
|
---|
| 228 | double sminl = 0;
|
---|
| 229 | double sminlo = 0;
|
---|
| 230 | double sminoa = 0;
|
---|
| 231 | double sn = 0;
|
---|
| 232 | double thresh = 0;
|
---|
| 233 | double tol = 0;
|
---|
| 234 | double tolmul = 0;
|
---|
| 235 | double unfl = 0;
|
---|
| 236 | double[] work0 = new double[0];
|
---|
| 237 | double[] work1 = new double[0];
|
---|
| 238 | double[] work2 = new double[0];
|
---|
| 239 | double[] work3 = new double[0];
|
---|
| 240 | int maxitr = 0;
|
---|
| 241 | bool matrixsplitflag = new bool();
|
---|
| 242 | bool iterflag = new bool();
|
---|
| 243 | double[] utemp = new double[0];
|
---|
| 244 | double[] vttemp = new double[0];
|
---|
| 245 | double[] ctemp = new double[0];
|
---|
| 246 | double[] etemp = new double[0];
|
---|
| 247 | bool rightside = new bool();
|
---|
| 248 | bool fwddir = new bool();
|
---|
| 249 | double tmp = 0;
|
---|
| 250 | int mm1 = 0;
|
---|
| 251 | int mm0 = 0;
|
---|
| 252 | bool bchangedir = new bool();
|
---|
| 253 | int uend = 0;
|
---|
| 254 | int cend = 0;
|
---|
| 255 | int vend = 0;
|
---|
| 256 | int i_ = 0;
|
---|
[2154] | 257 |
|
---|
[2430] | 258 | e = (double[])e.Clone();
|
---|
[2154] | 259 |
|
---|
[2430] | 260 | result = true;
|
---|
| 261 | if( n==0 )
|
---|
[2154] | 262 | {
|
---|
[2430] | 263 | return result;
|
---|
| 264 | }
|
---|
| 265 | if( n==1 )
|
---|
| 266 | {
|
---|
[2563] | 267 | if( (double)(d[1])<(double)(0) )
|
---|
[2154] | 268 | {
|
---|
[2430] | 269 | d[1] = -d[1];
|
---|
| 270 | if( ncvt>0 )
|
---|
[2154] | 271 | {
|
---|
[2430] | 272 | for(i_=vstart; i_<=vstart+ncvt-1;i_++)
|
---|
| 273 | {
|
---|
| 274 | vt[vstart,i_] = -1*vt[vstart,i_];
|
---|
| 275 | }
|
---|
[2154] | 276 | }
|
---|
| 277 | }
|
---|
[2430] | 278 | return result;
|
---|
[2154] | 279 | }
|
---|
| 280 |
|
---|
| 281 | //
|
---|
[2430] | 282 | // init
|
---|
[2154] | 283 | //
|
---|
[2430] | 284 | work0 = new double[n-1+1];
|
---|
| 285 | work1 = new double[n-1+1];
|
---|
| 286 | work2 = new double[n-1+1];
|
---|
| 287 | work3 = new double[n-1+1];
|
---|
| 288 | uend = ustart+Math.Max(nru-1, 0);
|
---|
| 289 | vend = vstart+Math.Max(ncvt-1, 0);
|
---|
| 290 | cend = cstart+Math.Max(ncc-1, 0);
|
---|
| 291 | utemp = new double[uend+1];
|
---|
| 292 | vttemp = new double[vend+1];
|
---|
| 293 | ctemp = new double[cend+1];
|
---|
| 294 | maxitr = 12;
|
---|
| 295 | rightside = true;
|
---|
| 296 | fwddir = true;
|
---|
| 297 |
|
---|
| 298 | //
|
---|
| 299 | // resize E from N-1 to N
|
---|
| 300 | //
|
---|
| 301 | etemp = new double[n+1];
|
---|
| 302 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 303 | {
|
---|
[2430] | 304 | etemp[i] = e[i];
|
---|
[2154] | 305 | }
|
---|
[2430] | 306 | e = new double[n+1];
|
---|
| 307 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 308 | {
|
---|
[2430] | 309 | e[i] = etemp[i];
|
---|
[2154] | 310 | }
|
---|
[2430] | 311 | e[n] = 0;
|
---|
| 312 | idir = 0;
|
---|
[2154] | 313 |
|
---|
| 314 | //
|
---|
[2430] | 315 | // Get machine constants
|
---|
[2154] | 316 | //
|
---|
[2430] | 317 | eps = AP.Math.MachineEpsilon;
|
---|
| 318 | unfl = AP.Math.MinRealNumber;
|
---|
| 319 |
|
---|
| 320 | //
|
---|
| 321 | // If matrix lower bidiagonal, rotate to be upper bidiagonal
|
---|
| 322 | // by applying Givens rotations on the left
|
---|
| 323 | //
|
---|
| 324 | if( !isupper )
|
---|
[2154] | 325 | {
|
---|
[2430] | 326 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 327 | {
|
---|
[2430] | 328 | rotations.generaterotation(d[i], e[i], ref cs, ref sn, ref r);
|
---|
| 329 | d[i] = r;
|
---|
| 330 | e[i] = sn*d[i+1];
|
---|
| 331 | d[i+1] = cs*d[i+1];
|
---|
| 332 | work0[i] = cs;
|
---|
| 333 | work1[i] = sn;
|
---|
[2154] | 334 | }
|
---|
[2430] | 335 |
|
---|
| 336 | //
|
---|
| 337 | // Update singular vectors if desired
|
---|
| 338 | //
|
---|
| 339 | if( nru>0 )
|
---|
| 340 | {
|
---|
| 341 | rotations.applyrotationsfromtheright(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, ref work0, ref work1, ref u, ref utemp);
|
---|
| 342 | }
|
---|
| 343 | if( ncc>0 )
|
---|
| 344 | {
|
---|
| 345 | rotations.applyrotationsfromtheleft(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
|
---|
| 346 | }
|
---|
[2154] | 347 | }
|
---|
| 348 |
|
---|
| 349 | //
|
---|
[2430] | 350 | // Compute singular values to relative accuracy TOL
|
---|
| 351 | // (By setting TOL to be negative, algorithm will compute
|
---|
| 352 | // singular values to absolute accuracy ABS(TOL)*norm(input matrix))
|
---|
[2154] | 353 | //
|
---|
[2430] | 354 | tolmul = Math.Max(10, Math.Min(100, Math.Pow(eps, -0.125)));
|
---|
| 355 | tol = tolmul*eps;
|
---|
| 356 | if( !isfractionalaccuracyrequired )
|
---|
| 357 | {
|
---|
| 358 | tol = -tol;
|
---|
| 359 | }
|
---|
[2154] | 360 |
|
---|
| 361 | //
|
---|
[2430] | 362 | // Compute approximate maximum, minimum singular values
|
---|
[2154] | 363 | //
|
---|
[2430] | 364 | smax = 0;
|
---|
| 365 | for(i=1; i<=n; i++)
|
---|
[2154] | 366 | {
|
---|
[2430] | 367 | smax = Math.Max(smax, Math.Abs(d[i]));
|
---|
[2154] | 368 | }
|
---|
[2430] | 369 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 370 | {
|
---|
[2430] | 371 | smax = Math.Max(smax, Math.Abs(e[i]));
|
---|
[2154] | 372 | }
|
---|
[2430] | 373 | sminl = 0;
|
---|
[2563] | 374 | if( (double)(tol)>=(double)(0) )
|
---|
[2154] | 375 | {
|
---|
[2430] | 376 |
|
---|
| 377 | //
|
---|
| 378 | // Relative accuracy desired
|
---|
| 379 | //
|
---|
| 380 | sminoa = Math.Abs(d[1]);
|
---|
[2563] | 381 | if( (double)(sminoa)!=(double)(0) )
|
---|
[2154] | 382 | {
|
---|
[2430] | 383 | mu = sminoa;
|
---|
| 384 | for(i=2; i<=n; i++)
|
---|
| 385 | {
|
---|
| 386 | mu = Math.Abs(d[i])*(mu/(mu+Math.Abs(e[i-1])));
|
---|
| 387 | sminoa = Math.Min(sminoa, mu);
|
---|
[2563] | 388 | if( (double)(sminoa)==(double)(0) )
|
---|
[2430] | 389 | {
|
---|
| 390 | break;
|
---|
| 391 | }
|
---|
| 392 | }
|
---|
[2154] | 393 | }
|
---|
[2430] | 394 | sminoa = sminoa/Math.Sqrt(n);
|
---|
| 395 | thresh = Math.Max(tol*sminoa, maxitr*n*n*unfl);
|
---|
[2154] | 396 | }
|
---|
| 397 | else
|
---|
| 398 | {
|
---|
| 399 |
|
---|
| 400 | //
|
---|
[2430] | 401 | // Absolute accuracy desired
|
---|
[2154] | 402 | //
|
---|
[2430] | 403 | thresh = Math.Max(Math.Abs(tol)*smax, maxitr*n*n*unfl);
|
---|
[2154] | 404 | }
|
---|
| 405 |
|
---|
| 406 | //
|
---|
[2430] | 407 | // Prepare for main iteration loop for the singular values
|
---|
| 408 | // (MAXIT is the maximum number of passes through the inner
|
---|
| 409 | // loop permitted before nonconvergence signalled.)
|
---|
[2154] | 410 | //
|
---|
[2430] | 411 | maxit = maxitr*n*n;
|
---|
| 412 | iter = 0;
|
---|
| 413 | oldll = -1;
|
---|
| 414 | oldm = -1;
|
---|
| 415 |
|
---|
| 416 | //
|
---|
| 417 | // M points to last element of unconverged part of matrix
|
---|
| 418 | //
|
---|
| 419 | m = n;
|
---|
| 420 |
|
---|
| 421 | //
|
---|
| 422 | // Begin main iteration loop
|
---|
| 423 | //
|
---|
| 424 | while( true )
|
---|
[2154] | 425 | {
|
---|
| 426 |
|
---|
| 427 | //
|
---|
[2430] | 428 | // Check for convergence or exceeding iteration count
|
---|
[2154] | 429 | //
|
---|
[2430] | 430 | if( m<=1 )
|
---|
| 431 | {
|
---|
| 432 | break;
|
---|
| 433 | }
|
---|
| 434 | if( iter>maxit )
|
---|
| 435 | {
|
---|
| 436 | result = false;
|
---|
| 437 | return result;
|
---|
| 438 | }
|
---|
[2154] | 439 |
|
---|
| 440 | //
|
---|
[2430] | 441 | // Find diagonal block of matrix to work on
|
---|
[2154] | 442 | //
|
---|
[2563] | 443 | if( (double)(tol)<(double)(0) & (double)(Math.Abs(d[m]))<=(double)(thresh) )
|
---|
[2154] | 444 | {
|
---|
[2430] | 445 | d[m] = 0;
|
---|
[2154] | 446 | }
|
---|
[2430] | 447 | smax = Math.Abs(d[m]);
|
---|
| 448 | smin = smax;
|
---|
| 449 | matrixsplitflag = false;
|
---|
| 450 | for(lll=1; lll<=m-1; lll++)
|
---|
[2154] | 451 | {
|
---|
[2430] | 452 | ll = m-lll;
|
---|
| 453 | abss = Math.Abs(d[ll]);
|
---|
| 454 | abse = Math.Abs(e[ll]);
|
---|
[2563] | 455 | if( (double)(tol)<(double)(0) & (double)(abss)<=(double)(thresh) )
|
---|
[2154] | 456 | {
|
---|
[2430] | 457 | d[ll] = 0;
|
---|
[2154] | 458 | }
|
---|
[2563] | 459 | if( (double)(abse)<=(double)(thresh) )
|
---|
[2154] | 460 | {
|
---|
[2430] | 461 | matrixsplitflag = true;
|
---|
| 462 | break;
|
---|
[2154] | 463 | }
|
---|
[2430] | 464 | smin = Math.Min(smin, abss);
|
---|
| 465 | smax = Math.Max(smax, Math.Max(abss, abse));
|
---|
[2154] | 466 | }
|
---|
[2430] | 467 | if( !matrixsplitflag )
|
---|
[2154] | 468 | {
|
---|
[2430] | 469 | ll = 0;
|
---|
[2154] | 470 | }
|
---|
| 471 | else
|
---|
| 472 | {
|
---|
| 473 |
|
---|
| 474 | //
|
---|
[2430] | 475 | // Matrix splits since E(LL) = 0
|
---|
[2154] | 476 | //
|
---|
[2430] | 477 | e[ll] = 0;
|
---|
| 478 | if( ll==m-1 )
|
---|
| 479 | {
|
---|
| 480 |
|
---|
| 481 | //
|
---|
| 482 | // Convergence of bottom singular value, return to top of loop
|
---|
| 483 | //
|
---|
| 484 | m = m-1;
|
---|
| 485 | continue;
|
---|
| 486 | }
|
---|
[2154] | 487 | }
|
---|
[2430] | 488 | ll = ll+1;
|
---|
[2154] | 489 |
|
---|
| 490 | //
|
---|
[2430] | 491 | // E(LL) through E(M-1) are nonzero, E(LL-1) is zero
|
---|
[2154] | 492 | //
|
---|
[2430] | 493 | if( ll==m-1 )
|
---|
[2154] | 494 | {
|
---|
[2430] | 495 |
|
---|
| 496 | //
|
---|
| 497 | // 2 by 2 block, handle separately
|
---|
| 498 | //
|
---|
| 499 | svdv2x2(d[m-1], e[m-1], d[m], ref sigmn, ref sigmx, ref sinr, ref cosr, ref sinl, ref cosl);
|
---|
| 500 | d[m-1] = sigmx;
|
---|
[2154] | 501 | e[m-1] = 0;
|
---|
[2430] | 502 | d[m] = sigmn;
|
---|
[2154] | 503 |
|
---|
| 504 | //
|
---|
[2430] | 505 | // Compute singular vectors, if desired
|
---|
[2154] | 506 | //
|
---|
[2430] | 507 | if( ncvt>0 )
|
---|
[2154] | 508 | {
|
---|
[2430] | 509 | mm0 = m+(vstart-1);
|
---|
| 510 | mm1 = m-1+(vstart-1);
|
---|
| 511 | for(i_=vstart; i_<=vend;i_++)
|
---|
[2154] | 512 | {
|
---|
[2430] | 513 | vttemp[i_] = cosr*vt[mm1,i_];
|
---|
[2154] | 514 | }
|
---|
[2430] | 515 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 516 | {
|
---|
| 517 | vttemp[i_] = vttemp[i_] + sinr*vt[mm0,i_];
|
---|
| 518 | }
|
---|
| 519 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 520 | {
|
---|
| 521 | vt[mm0,i_] = cosr*vt[mm0,i_];
|
---|
| 522 | }
|
---|
| 523 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 524 | {
|
---|
| 525 | vt[mm0,i_] = vt[mm0,i_] - sinr*vt[mm1,i_];
|
---|
| 526 | }
|
---|
| 527 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 528 | {
|
---|
| 529 | vt[mm1,i_] = vttemp[i_];
|
---|
| 530 | }
|
---|
[2154] | 531 | }
|
---|
[2430] | 532 | if( nru>0 )
|
---|
[2154] | 533 | {
|
---|
[2430] | 534 | mm0 = m+ustart-1;
|
---|
| 535 | mm1 = m-1+ustart-1;
|
---|
| 536 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 537 | {
|
---|
| 538 | utemp[i_] = cosl*u[i_,mm1];
|
---|
| 539 | }
|
---|
| 540 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 541 | {
|
---|
| 542 | utemp[i_] = utemp[i_] + sinl*u[i_,mm0];
|
---|
| 543 | }
|
---|
| 544 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 545 | {
|
---|
| 546 | u[i_,mm0] = cosl*u[i_,mm0];
|
---|
| 547 | }
|
---|
| 548 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 549 | {
|
---|
| 550 | u[i_,mm0] = u[i_,mm0] - sinl*u[i_,mm1];
|
---|
| 551 | }
|
---|
| 552 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 553 | {
|
---|
| 554 | u[i_,mm1] = utemp[i_];
|
---|
| 555 | }
|
---|
[2154] | 556 | }
|
---|
[2430] | 557 | if( ncc>0 )
|
---|
[2154] | 558 | {
|
---|
[2430] | 559 | mm0 = m+cstart-1;
|
---|
| 560 | mm1 = m-1+cstart-1;
|
---|
| 561 | for(i_=cstart; i_<=cend;i_++)
|
---|
[2154] | 562 | {
|
---|
[2430] | 563 | ctemp[i_] = cosl*c[mm1,i_];
|
---|
[2154] | 564 | }
|
---|
[2430] | 565 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 566 | {
|
---|
| 567 | ctemp[i_] = ctemp[i_] + sinl*c[mm0,i_];
|
---|
| 568 | }
|
---|
| 569 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 570 | {
|
---|
| 571 | c[mm0,i_] = cosl*c[mm0,i_];
|
---|
| 572 | }
|
---|
| 573 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 574 | {
|
---|
| 575 | c[mm0,i_] = c[mm0,i_] - sinl*c[mm1,i_];
|
---|
| 576 | }
|
---|
| 577 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 578 | {
|
---|
| 579 | c[mm1,i_] = ctemp[i_];
|
---|
| 580 | }
|
---|
[2154] | 581 | }
|
---|
[2430] | 582 | m = m-2;
|
---|
| 583 | continue;
|
---|
[2154] | 584 | }
|
---|
| 585 |
|
---|
| 586 | //
|
---|
[2430] | 587 | // If working on new submatrix, choose shift direction
|
---|
| 588 | // (from larger end diagonal element towards smaller)
|
---|
[2154] | 589 | //
|
---|
[2430] | 590 | // Previously was
|
---|
| 591 | // "if (LL>OLDM) or (M<OLDLL) then"
|
---|
| 592 | // fixed thanks to Michael Rolle < m@rolle.name >
|
---|
| 593 | // Very strange that LAPACK still contains it.
|
---|
[2154] | 594 | //
|
---|
[2430] | 595 | bchangedir = false;
|
---|
[2563] | 596 | if( idir==1 & (double)(Math.Abs(d[ll]))<(double)(1.0E-3*Math.Abs(d[m])) )
|
---|
[2154] | 597 | {
|
---|
[2430] | 598 | bchangedir = true;
|
---|
[2154] | 599 | }
|
---|
[2563] | 600 | if( idir==2 & (double)(Math.Abs(d[m]))<(double)(1.0E-3*Math.Abs(d[ll])) )
|
---|
[2154] | 601 | {
|
---|
[2430] | 602 | bchangedir = true;
|
---|
[2154] | 603 | }
|
---|
[2430] | 604 | if( ll!=oldll | m!=oldm | bchangedir )
|
---|
[2154] | 605 | {
|
---|
[2563] | 606 | if( (double)(Math.Abs(d[ll]))>=(double)(Math.Abs(d[m])) )
|
---|
[2154] | 607 | {
|
---|
[2430] | 608 |
|
---|
| 609 | //
|
---|
| 610 | // Chase bulge from top (big end) to bottom (small end)
|
---|
| 611 | //
|
---|
| 612 | idir = 1;
|
---|
[2154] | 613 | }
|
---|
[2430] | 614 | else
|
---|
| 615 | {
|
---|
| 616 |
|
---|
| 617 | //
|
---|
| 618 | // Chase bulge from bottom (big end) to top (small end)
|
---|
| 619 | //
|
---|
| 620 | idir = 2;
|
---|
| 621 | }
|
---|
[2154] | 622 | }
|
---|
[2430] | 623 |
|
---|
| 624 | //
|
---|
| 625 | // Apply convergence tests
|
---|
| 626 | //
|
---|
[2154] | 627 | if( idir==1 )
|
---|
| 628 | {
|
---|
| 629 |
|
---|
| 630 | //
|
---|
[2430] | 631 | // Run convergence test in forward direction
|
---|
| 632 | // First apply standard test to bottom of matrix
|
---|
[2154] | 633 | //
|
---|
[2563] | 634 | if( (double)(Math.Abs(e[m-1]))<=(double)(Math.Abs(tol)*Math.Abs(d[m])) | (double)(tol)<(double)(0) & (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
|
---|
[2154] | 635 | {
|
---|
[2430] | 636 | e[m-1] = 0;
|
---|
| 637 | continue;
|
---|
| 638 | }
|
---|
[2563] | 639 | if( (double)(tol)>=(double)(0) )
|
---|
[2430] | 640 | {
|
---|
| 641 |
|
---|
| 642 | //
|
---|
| 643 | // If relative accuracy desired,
|
---|
| 644 | // apply convergence criterion forward
|
---|
| 645 | //
|
---|
| 646 | mu = Math.Abs(d[ll]);
|
---|
| 647 | sminl = mu;
|
---|
| 648 | iterflag = false;
|
---|
| 649 | for(lll=ll; lll<=m-1; lll++)
|
---|
[2154] | 650 | {
|
---|
[2563] | 651 | if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
|
---|
[2430] | 652 | {
|
---|
| 653 | e[lll] = 0;
|
---|
| 654 | iterflag = true;
|
---|
| 655 | break;
|
---|
| 656 | }
|
---|
| 657 | sminlo = sminl;
|
---|
| 658 | mu = Math.Abs(d[lll+1])*(mu/(mu+Math.Abs(e[lll])));
|
---|
| 659 | sminl = Math.Min(sminl, mu);
|
---|
[2154] | 660 | }
|
---|
[2430] | 661 | if( iterflag )
|
---|
| 662 | {
|
---|
| 663 | continue;
|
---|
| 664 | }
|
---|
[2154] | 665 | }
|
---|
[2430] | 666 | }
|
---|
| 667 | else
|
---|
| 668 | {
|
---|
[2154] | 669 |
|
---|
| 670 | //
|
---|
[2430] | 671 | // Run convergence test in backward direction
|
---|
| 672 | // First apply standard test to top of matrix
|
---|
[2154] | 673 | //
|
---|
[2563] | 674 | if( (double)(Math.Abs(e[ll]))<=(double)(Math.Abs(tol)*Math.Abs(d[ll])) | (double)(tol)<(double)(0) & (double)(Math.Abs(e[ll]))<=(double)(thresh) )
|
---|
[2154] | 675 | {
|
---|
[2430] | 676 | e[ll] = 0;
|
---|
| 677 | continue;
|
---|
[2154] | 678 | }
|
---|
[2563] | 679 | if( (double)(tol)>=(double)(0) )
|
---|
[2154] | 680 | {
|
---|
[2430] | 681 |
|
---|
| 682 | //
|
---|
| 683 | // If relative accuracy desired,
|
---|
| 684 | // apply convergence criterion backward
|
---|
| 685 | //
|
---|
| 686 | mu = Math.Abs(d[m]);
|
---|
| 687 | sminl = mu;
|
---|
| 688 | iterflag = false;
|
---|
| 689 | for(lll=m-1; lll>=ll; lll--)
|
---|
| 690 | {
|
---|
[2563] | 691 | if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
|
---|
[2430] | 692 | {
|
---|
| 693 | e[lll] = 0;
|
---|
| 694 | iterflag = true;
|
---|
| 695 | break;
|
---|
| 696 | }
|
---|
| 697 | sminlo = sminl;
|
---|
| 698 | mu = Math.Abs(d[lll])*(mu/(mu+Math.Abs(e[lll])));
|
---|
| 699 | sminl = Math.Min(sminl, mu);
|
---|
| 700 | }
|
---|
| 701 | if( iterflag )
|
---|
| 702 | {
|
---|
| 703 | continue;
|
---|
| 704 | }
|
---|
[2154] | 705 | }
|
---|
[2430] | 706 | }
|
---|
| 707 | oldll = ll;
|
---|
| 708 | oldm = m;
|
---|
| 709 |
|
---|
| 710 | //
|
---|
| 711 | // Compute shift. First, test if shifting would ruin relative
|
---|
| 712 | // accuracy, and if so set the shift to zero.
|
---|
| 713 | //
|
---|
[2563] | 714 | if( (double)(tol)>=(double)(0) & (double)(n*tol*(sminl/smax))<=(double)(Math.Max(eps, 0.01*tol)) )
|
---|
[2430] | 715 | {
|
---|
[2154] | 716 |
|
---|
| 717 | //
|
---|
[2430] | 718 | // Use a zero shift to avoid loss of relative accuracy
|
---|
[2154] | 719 | //
|
---|
[2430] | 720 | shift = 0;
|
---|
[2154] | 721 | }
|
---|
| 722 | else
|
---|
| 723 | {
|
---|
| 724 |
|
---|
| 725 | //
|
---|
[2430] | 726 | // Compute the shift from 2-by-2 block at end of matrix
|
---|
[2154] | 727 | //
|
---|
[2430] | 728 | if( idir==1 )
|
---|
[2154] | 729 | {
|
---|
[2430] | 730 | sll = Math.Abs(d[ll]);
|
---|
| 731 | svd2x2(d[m-1], e[m-1], d[m], ref shift, ref r);
|
---|
[2154] | 732 | }
|
---|
[2430] | 733 | else
|
---|
[2154] | 734 | {
|
---|
[2430] | 735 | sll = Math.Abs(d[m]);
|
---|
| 736 | svd2x2(d[ll], e[ll], d[ll+1], ref shift, ref r);
|
---|
[2154] | 737 | }
|
---|
| 738 |
|
---|
| 739 | //
|
---|
[2430] | 740 | // Test if shift negligible, and if so set to zero
|
---|
[2154] | 741 | //
|
---|
[2563] | 742 | if( (double)(sll)>(double)(0) )
|
---|
[2154] | 743 | {
|
---|
[2563] | 744 | if( (double)(AP.Math.Sqr(shift/sll))<(double)(eps) )
|
---|
[2430] | 745 | {
|
---|
| 746 | shift = 0;
|
---|
| 747 | }
|
---|
[2154] | 748 | }
|
---|
| 749 | }
|
---|
| 750 |
|
---|
| 751 | //
|
---|
[2430] | 752 | // Increment iteration count
|
---|
[2154] | 753 | //
|
---|
[2430] | 754 | iter = iter+m-ll;
|
---|
| 755 |
|
---|
| 756 | //
|
---|
| 757 | // If SHIFT = 0, do simplified QR iteration
|
---|
| 758 | //
|
---|
[2563] | 759 | if( (double)(shift)==(double)(0) )
|
---|
[2154] | 760 | {
|
---|
[2430] | 761 | if( idir==1 )
|
---|
[2154] | 762 | {
|
---|
[2430] | 763 |
|
---|
| 764 | //
|
---|
| 765 | // Chase bulge from top to bottom
|
---|
| 766 | // Save cosines and sines for later singular vector updates
|
---|
| 767 | //
|
---|
| 768 | cs = 1;
|
---|
| 769 | oldcs = 1;
|
---|
| 770 | for(i=ll; i<=m-1; i++)
|
---|
[2154] | 771 | {
|
---|
[2430] | 772 | rotations.generaterotation(d[i]*cs, e[i], ref cs, ref sn, ref r);
|
---|
| 773 | if( i>ll )
|
---|
| 774 | {
|
---|
| 775 | e[i-1] = oldsn*r;
|
---|
| 776 | }
|
---|
| 777 | rotations.generaterotation(oldcs*r, d[i+1]*sn, ref oldcs, ref oldsn, ref tmp);
|
---|
| 778 | d[i] = tmp;
|
---|
| 779 | work0[i-ll+1] = cs;
|
---|
| 780 | work1[i-ll+1] = sn;
|
---|
| 781 | work2[i-ll+1] = oldcs;
|
---|
| 782 | work3[i-ll+1] = oldsn;
|
---|
[2154] | 783 | }
|
---|
[2430] | 784 | h = d[m]*cs;
|
---|
| 785 | d[m] = h*oldcs;
|
---|
| 786 | e[m-1] = h*oldsn;
|
---|
| 787 |
|
---|
| 788 | //
|
---|
| 789 | // Update singular vectors
|
---|
| 790 | //
|
---|
| 791 | if( ncvt>0 )
|
---|
[2154] | 792 | {
|
---|
[2430] | 793 | rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
|
---|
[2154] | 794 | }
|
---|
[2430] | 795 | if( nru>0 )
|
---|
| 796 | {
|
---|
| 797 | rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
|
---|
| 798 | }
|
---|
| 799 | if( ncc>0 )
|
---|
| 800 | {
|
---|
| 801 | rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
|
---|
| 802 | }
|
---|
| 803 |
|
---|
| 804 | //
|
---|
| 805 | // Test convergence
|
---|
| 806 | //
|
---|
[2563] | 807 | if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
|
---|
[2430] | 808 | {
|
---|
| 809 | e[m-1] = 0;
|
---|
| 810 | }
|
---|
[2154] | 811 | }
|
---|
[2430] | 812 | else
|
---|
[2154] | 813 | {
|
---|
[2430] | 814 |
|
---|
| 815 | //
|
---|
| 816 | // Chase bulge from bottom to top
|
---|
| 817 | // Save cosines and sines for later singular vector updates
|
---|
| 818 | //
|
---|
| 819 | cs = 1;
|
---|
| 820 | oldcs = 1;
|
---|
| 821 | for(i=m; i>=ll+1; i--)
|
---|
| 822 | {
|
---|
| 823 | rotations.generaterotation(d[i]*cs, e[i-1], ref cs, ref sn, ref r);
|
---|
| 824 | if( i<m )
|
---|
| 825 | {
|
---|
| 826 | e[i] = oldsn*r;
|
---|
| 827 | }
|
---|
| 828 | rotations.generaterotation(oldcs*r, d[i-1]*sn, ref oldcs, ref oldsn, ref tmp);
|
---|
| 829 | d[i] = tmp;
|
---|
| 830 | work0[i-ll] = cs;
|
---|
| 831 | work1[i-ll] = -sn;
|
---|
| 832 | work2[i-ll] = oldcs;
|
---|
| 833 | work3[i-ll] = -oldsn;
|
---|
| 834 | }
|
---|
| 835 | h = d[ll]*cs;
|
---|
| 836 | d[ll] = h*oldcs;
|
---|
| 837 | e[ll] = h*oldsn;
|
---|
| 838 |
|
---|
| 839 | //
|
---|
| 840 | // Update singular vectors
|
---|
| 841 | //
|
---|
| 842 | if( ncvt>0 )
|
---|
| 843 | {
|
---|
| 844 | rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
|
---|
| 845 | }
|
---|
| 846 | if( nru>0 )
|
---|
| 847 | {
|
---|
| 848 | rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
|
---|
| 849 | }
|
---|
| 850 | if( ncc>0 )
|
---|
| 851 | {
|
---|
| 852 | rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
|
---|
| 853 | }
|
---|
| 854 |
|
---|
| 855 | //
|
---|
| 856 | // Test convergence
|
---|
| 857 | //
|
---|
[2563] | 858 | if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
|
---|
[2430] | 859 | {
|
---|
| 860 | e[ll] = 0;
|
---|
| 861 | }
|
---|
[2154] | 862 | }
|
---|
| 863 | }
|
---|
| 864 | else
|
---|
| 865 | {
|
---|
| 866 |
|
---|
| 867 | //
|
---|
[2430] | 868 | // Use nonzero shift
|
---|
[2154] | 869 | //
|
---|
[2430] | 870 | if( idir==1 )
|
---|
[2154] | 871 | {
|
---|
[2430] | 872 |
|
---|
| 873 | //
|
---|
| 874 | // Chase bulge from top to bottom
|
---|
| 875 | // Save cosines and sines for later singular vector updates
|
---|
| 876 | //
|
---|
| 877 | f = (Math.Abs(d[ll])-shift)*(extsignbdsqr(1, d[ll])+shift/d[ll]);
|
---|
| 878 | g = e[ll];
|
---|
| 879 | for(i=ll; i<=m-1; i++)
|
---|
[2154] | 880 | {
|
---|
[2430] | 881 | rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
|
---|
| 882 | if( i>ll )
|
---|
| 883 | {
|
---|
| 884 | e[i-1] = r;
|
---|
| 885 | }
|
---|
| 886 | f = cosr*d[i]+sinr*e[i];
|
---|
| 887 | e[i] = cosr*e[i]-sinr*d[i];
|
---|
| 888 | g = sinr*d[i+1];
|
---|
| 889 | d[i+1] = cosr*d[i+1];
|
---|
| 890 | rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
|
---|
| 891 | d[i] = r;
|
---|
| 892 | f = cosl*e[i]+sinl*d[i+1];
|
---|
| 893 | d[i+1] = cosl*d[i+1]-sinl*e[i];
|
---|
| 894 | if( i<m-1 )
|
---|
| 895 | {
|
---|
| 896 | g = sinl*e[i+1];
|
---|
| 897 | e[i+1] = cosl*e[i+1];
|
---|
| 898 | }
|
---|
| 899 | work0[i-ll+1] = cosr;
|
---|
| 900 | work1[i-ll+1] = sinr;
|
---|
| 901 | work2[i-ll+1] = cosl;
|
---|
| 902 | work3[i-ll+1] = sinl;
|
---|
[2154] | 903 | }
|
---|
[2430] | 904 | e[m-1] = f;
|
---|
| 905 |
|
---|
| 906 | //
|
---|
| 907 | // Update singular vectors
|
---|
| 908 | //
|
---|
| 909 | if( ncvt>0 )
|
---|
[2154] | 910 | {
|
---|
[2430] | 911 | rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
|
---|
[2154] | 912 | }
|
---|
[2430] | 913 | if( nru>0 )
|
---|
| 914 | {
|
---|
| 915 | rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
|
---|
| 916 | }
|
---|
| 917 | if( ncc>0 )
|
---|
| 918 | {
|
---|
| 919 | rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
|
---|
| 920 | }
|
---|
| 921 |
|
---|
| 922 | //
|
---|
| 923 | // Test convergence
|
---|
| 924 | //
|
---|
[2563] | 925 | if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
|
---|
[2430] | 926 | {
|
---|
| 927 | e[m-1] = 0;
|
---|
| 928 | }
|
---|
[2154] | 929 | }
|
---|
[2430] | 930 | else
|
---|
[2154] | 931 | {
|
---|
[2430] | 932 |
|
---|
| 933 | //
|
---|
| 934 | // Chase bulge from bottom to top
|
---|
| 935 | // Save cosines and sines for later singular vector updates
|
---|
| 936 | //
|
---|
| 937 | f = (Math.Abs(d[m])-shift)*(extsignbdsqr(1, d[m])+shift/d[m]);
|
---|
| 938 | g = e[m-1];
|
---|
| 939 | for(i=m; i>=ll+1; i--)
|
---|
| 940 | {
|
---|
| 941 | rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
|
---|
| 942 | if( i<m )
|
---|
| 943 | {
|
---|
| 944 | e[i] = r;
|
---|
| 945 | }
|
---|
| 946 | f = cosr*d[i]+sinr*e[i-1];
|
---|
| 947 | e[i-1] = cosr*e[i-1]-sinr*d[i];
|
---|
| 948 | g = sinr*d[i-1];
|
---|
| 949 | d[i-1] = cosr*d[i-1];
|
---|
| 950 | rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
|
---|
| 951 | d[i] = r;
|
---|
| 952 | f = cosl*e[i-1]+sinl*d[i-1];
|
---|
| 953 | d[i-1] = cosl*d[i-1]-sinl*e[i-1];
|
---|
| 954 | if( i>ll+1 )
|
---|
| 955 | {
|
---|
| 956 | g = sinl*e[i-2];
|
---|
| 957 | e[i-2] = cosl*e[i-2];
|
---|
| 958 | }
|
---|
| 959 | work0[i-ll] = cosr;
|
---|
| 960 | work1[i-ll] = -sinr;
|
---|
| 961 | work2[i-ll] = cosl;
|
---|
| 962 | work3[i-ll] = -sinl;
|
---|
| 963 | }
|
---|
| 964 | e[ll] = f;
|
---|
| 965 |
|
---|
| 966 | //
|
---|
| 967 | // Test convergence
|
---|
| 968 | //
|
---|
[2563] | 969 | if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
|
---|
[2430] | 970 | {
|
---|
| 971 | e[ll] = 0;
|
---|
| 972 | }
|
---|
| 973 |
|
---|
| 974 | //
|
---|
| 975 | // Update singular vectors if desired
|
---|
| 976 | //
|
---|
| 977 | if( ncvt>0 )
|
---|
| 978 | {
|
---|
| 979 | rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
|
---|
| 980 | }
|
---|
| 981 | if( nru>0 )
|
---|
| 982 | {
|
---|
| 983 | rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
|
---|
| 984 | }
|
---|
| 985 | if( ncc>0 )
|
---|
| 986 | {
|
---|
| 987 | rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
|
---|
| 988 | }
|
---|
[2154] | 989 | }
|
---|
| 990 | }
|
---|
[2430] | 991 |
|
---|
| 992 | //
|
---|
| 993 | // QR iteration finished, go back and check convergence
|
---|
| 994 | //
|
---|
| 995 | continue;
|
---|
[2154] | 996 | }
|
---|
| 997 |
|
---|
| 998 | //
|
---|
[2430] | 999 | // All singular values converged, so make them positive
|
---|
[2154] | 1000 | //
|
---|
[2430] | 1001 | for(i=1; i<=n; i++)
|
---|
[2154] | 1002 | {
|
---|
[2563] | 1003 | if( (double)(d[i])<(double)(0) )
|
---|
[2154] | 1004 | {
|
---|
[2430] | 1005 | d[i] = -d[i];
|
---|
| 1006 |
|
---|
| 1007 | //
|
---|
| 1008 | // Change sign of singular vectors, if desired
|
---|
| 1009 | //
|
---|
| 1010 | if( ncvt>0 )
|
---|
[2154] | 1011 | {
|
---|
[2430] | 1012 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 1013 | {
|
---|
| 1014 | vt[i+vstart-1,i_] = -1*vt[i+vstart-1,i_];
|
---|
| 1015 | }
|
---|
[2154] | 1016 | }
|
---|
| 1017 | }
|
---|
| 1018 | }
|
---|
| 1019 |
|
---|
| 1020 | //
|
---|
[2430] | 1021 | // Sort the singular values into decreasing order (insertion sort on
|
---|
| 1022 | // singular values, but only one transposition per singular vector)
|
---|
[2154] | 1023 | //
|
---|
[2430] | 1024 | for(i=1; i<=n-1; i++)
|
---|
[2154] | 1025 | {
|
---|
| 1026 |
|
---|
| 1027 | //
|
---|
[2430] | 1028 | // Scan for smallest D(I)
|
---|
[2154] | 1029 | //
|
---|
[2430] | 1030 | isub = 1;
|
---|
| 1031 | smin = d[1];
|
---|
| 1032 | for(j=2; j<=n+1-i; j++)
|
---|
[2154] | 1033 | {
|
---|
[2563] | 1034 | if( (double)(d[j])<=(double)(smin) )
|
---|
[2154] | 1035 | {
|
---|
[2430] | 1036 | isub = j;
|
---|
| 1037 | smin = d[j];
|
---|
[2154] | 1038 | }
|
---|
| 1039 | }
|
---|
[2430] | 1040 | if( isub!=n+1-i )
|
---|
[2154] | 1041 | {
|
---|
[2430] | 1042 |
|
---|
| 1043 | //
|
---|
| 1044 | // Swap singular values and vectors
|
---|
| 1045 | //
|
---|
| 1046 | d[isub] = d[n+1-i];
|
---|
| 1047 | d[n+1-i] = smin;
|
---|
| 1048 | if( ncvt>0 )
|
---|
[2154] | 1049 | {
|
---|
[2430] | 1050 | j = n+1-i;
|
---|
| 1051 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 1052 | {
|
---|
| 1053 | vttemp[i_] = vt[isub+vstart-1,i_];
|
---|
| 1054 | }
|
---|
| 1055 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 1056 | {
|
---|
| 1057 | vt[isub+vstart-1,i_] = vt[j+vstart-1,i_];
|
---|
| 1058 | }
|
---|
| 1059 | for(i_=vstart; i_<=vend;i_++)
|
---|
| 1060 | {
|
---|
| 1061 | vt[j+vstart-1,i_] = vttemp[i_];
|
---|
| 1062 | }
|
---|
[2154] | 1063 | }
|
---|
[2430] | 1064 | if( nru>0 )
|
---|
[2154] | 1065 | {
|
---|
[2430] | 1066 | j = n+1-i;
|
---|
| 1067 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 1068 | {
|
---|
| 1069 | utemp[i_] = u[i_,isub+ustart-1];
|
---|
| 1070 | }
|
---|
| 1071 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 1072 | {
|
---|
| 1073 | u[i_,isub+ustart-1] = u[i_,j+ustart-1];
|
---|
| 1074 | }
|
---|
| 1075 | for(i_=ustart; i_<=uend;i_++)
|
---|
| 1076 | {
|
---|
| 1077 | u[i_,j+ustart-1] = utemp[i_];
|
---|
| 1078 | }
|
---|
[2154] | 1079 | }
|
---|
[2430] | 1080 | if( ncc>0 )
|
---|
[2154] | 1081 | {
|
---|
[2430] | 1082 | j = n+1-i;
|
---|
| 1083 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 1084 | {
|
---|
| 1085 | ctemp[i_] = c[isub+cstart-1,i_];
|
---|
| 1086 | }
|
---|
| 1087 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 1088 | {
|
---|
| 1089 | c[isub+cstart-1,i_] = c[j+cstart-1,i_];
|
---|
| 1090 | }
|
---|
| 1091 | for(i_=cstart; i_<=cend;i_++)
|
---|
| 1092 | {
|
---|
| 1093 | c[j+cstart-1,i_] = ctemp[i_];
|
---|
| 1094 | }
|
---|
[2154] | 1095 | }
|
---|
| 1096 | }
|
---|
| 1097 | }
|
---|
[2430] | 1098 | return result;
|
---|
[2154] | 1099 | }
|
---|
| 1100 |
|
---|
| 1101 |
|
---|
[2430] | 1102 | private static double extsignbdsqr(double a,
|
---|
| 1103 | double b)
|
---|
[2154] | 1104 | {
|
---|
[2430] | 1105 | double result = 0;
|
---|
[2154] | 1106 |
|
---|
[2563] | 1107 | if( (double)(b)>=(double)(0) )
|
---|
[2154] | 1108 | {
|
---|
[2430] | 1109 | result = Math.Abs(a);
|
---|
[2154] | 1110 | }
|
---|
| 1111 | else
|
---|
| 1112 | {
|
---|
[2430] | 1113 | result = -Math.Abs(a);
|
---|
[2154] | 1114 | }
|
---|
[2430] | 1115 | return result;
|
---|
[2154] | 1116 | }
|
---|
[2430] | 1117 |
|
---|
| 1118 |
|
---|
| 1119 | private static void svd2x2(double f,
|
---|
| 1120 | double g,
|
---|
| 1121 | double h,
|
---|
| 1122 | ref double ssmin,
|
---|
| 1123 | ref double ssmax)
|
---|
[2154] | 1124 | {
|
---|
[2430] | 1125 | double aas = 0;
|
---|
| 1126 | double at = 0;
|
---|
| 1127 | double au = 0;
|
---|
| 1128 | double c = 0;
|
---|
| 1129 | double fa = 0;
|
---|
| 1130 | double fhmn = 0;
|
---|
| 1131 | double fhmx = 0;
|
---|
| 1132 | double ga = 0;
|
---|
| 1133 | double ha = 0;
|
---|
| 1134 |
|
---|
| 1135 | fa = Math.Abs(f);
|
---|
| 1136 | ga = Math.Abs(g);
|
---|
| 1137 | ha = Math.Abs(h);
|
---|
| 1138 | fhmn = Math.Min(fa, ha);
|
---|
| 1139 | fhmx = Math.Max(fa, ha);
|
---|
[2563] | 1140 | if( (double)(fhmn)==(double)(0) )
|
---|
[2154] | 1141 | {
|
---|
[2430] | 1142 | ssmin = 0;
|
---|
[2563] | 1143 | if( (double)(fhmx)==(double)(0) )
|
---|
[2154] | 1144 | {
|
---|
| 1145 | ssmax = ga;
|
---|
| 1146 | }
|
---|
| 1147 | else
|
---|
| 1148 | {
|
---|
[2430] | 1149 | ssmax = Math.Max(fhmx, ga)*Math.Sqrt(1+AP.Math.Sqr(Math.Min(fhmx, ga)/Math.Max(fhmx, ga)));
|
---|
| 1150 | }
|
---|
| 1151 | }
|
---|
| 1152 | else
|
---|
| 1153 | {
|
---|
[2563] | 1154 | if( (double)(ga)<(double)(fhmx) )
|
---|
[2430] | 1155 | {
|
---|
[2154] | 1156 | aas = 1+fhmn/fhmx;
|
---|
| 1157 | at = (fhmx-fhmn)/fhmx;
|
---|
[2430] | 1158 | au = AP.Math.Sqr(ga/fhmx);
|
---|
| 1159 | c = 2/(Math.Sqrt(aas*aas+au)+Math.Sqrt(at*at+au));
|
---|
| 1160 | ssmin = fhmn*c;
|
---|
| 1161 | ssmax = fhmx/c;
|
---|
[2154] | 1162 | }
|
---|
[2430] | 1163 | else
|
---|
| 1164 | {
|
---|
| 1165 | au = fhmx/ga;
|
---|
[2563] | 1166 | if( (double)(au)==(double)(0) )
|
---|
[2430] | 1167 | {
|
---|
| 1168 |
|
---|
| 1169 | //
|
---|
| 1170 | // Avoid possible harmful underflow if exponent range
|
---|
| 1171 | // asymmetric (true SSMIN may not underflow even if
|
---|
| 1172 | // AU underflows)
|
---|
| 1173 | //
|
---|
| 1174 | ssmin = fhmn*fhmx/ga;
|
---|
| 1175 | ssmax = ga;
|
---|
| 1176 | }
|
---|
| 1177 | else
|
---|
| 1178 | {
|
---|
| 1179 | aas = 1+fhmn/fhmx;
|
---|
| 1180 | at = (fhmx-fhmn)/fhmx;
|
---|
| 1181 | c = 1/(Math.Sqrt(1+AP.Math.Sqr(aas*au))+Math.Sqrt(1+AP.Math.Sqr(at*au)));
|
---|
| 1182 | ssmin = fhmn*c*au;
|
---|
| 1183 | ssmin = ssmin+ssmin;
|
---|
| 1184 | ssmax = ga/(c+c);
|
---|
| 1185 | }
|
---|
| 1186 | }
|
---|
[2154] | 1187 | }
|
---|
| 1188 | }
|
---|
| 1189 |
|
---|
| 1190 |
|
---|
[2430] | 1191 | private static void svdv2x2(double f,
|
---|
| 1192 | double g,
|
---|
| 1193 | double h,
|
---|
| 1194 | ref double ssmin,
|
---|
| 1195 | ref double ssmax,
|
---|
| 1196 | ref double snr,
|
---|
| 1197 | ref double csr,
|
---|
| 1198 | ref double snl,
|
---|
| 1199 | ref double csl)
|
---|
| 1200 | {
|
---|
| 1201 | bool gasmal = new bool();
|
---|
| 1202 | bool swp = new bool();
|
---|
| 1203 | int pmax = 0;
|
---|
| 1204 | double a = 0;
|
---|
| 1205 | double clt = 0;
|
---|
| 1206 | double crt = 0;
|
---|
| 1207 | double d = 0;
|
---|
| 1208 | double fa = 0;
|
---|
| 1209 | double ft = 0;
|
---|
| 1210 | double ga = 0;
|
---|
| 1211 | double gt = 0;
|
---|
| 1212 | double ha = 0;
|
---|
| 1213 | double ht = 0;
|
---|
| 1214 | double l = 0;
|
---|
| 1215 | double m = 0;
|
---|
| 1216 | double mm = 0;
|
---|
| 1217 | double r = 0;
|
---|
| 1218 | double s = 0;
|
---|
| 1219 | double slt = 0;
|
---|
| 1220 | double srt = 0;
|
---|
| 1221 | double t = 0;
|
---|
| 1222 | double temp = 0;
|
---|
| 1223 | double tsign = 0;
|
---|
| 1224 | double tt = 0;
|
---|
| 1225 | double v = 0;
|
---|
[2154] | 1226 |
|
---|
[2430] | 1227 | ft = f;
|
---|
| 1228 | fa = Math.Abs(ft);
|
---|
| 1229 | ht = h;
|
---|
| 1230 | ha = Math.Abs(h);
|
---|
[2154] | 1231 |
|
---|
| 1232 | //
|
---|
[2430] | 1233 | // PMAX points to the maximum absolute element of matrix
|
---|
| 1234 | // PMAX = 1 if F largest in absolute values
|
---|
| 1235 | // PMAX = 2 if G largest in absolute values
|
---|
| 1236 | // PMAX = 3 if H largest in absolute values
|
---|
[2154] | 1237 | //
|
---|
[2430] | 1238 | pmax = 1;
|
---|
[2563] | 1239 | swp = (double)(ha)>(double)(fa);
|
---|
[2430] | 1240 | if( swp )
|
---|
[2154] | 1241 | {
|
---|
[2430] | 1242 |
|
---|
| 1243 | //
|
---|
| 1244 | // Now FA .ge. HA
|
---|
| 1245 | //
|
---|
| 1246 | pmax = 3;
|
---|
| 1247 | temp = ft;
|
---|
| 1248 | ft = ht;
|
---|
| 1249 | ht = temp;
|
---|
| 1250 | temp = fa;
|
---|
| 1251 | fa = ha;
|
---|
| 1252 | ha = temp;
|
---|
[2154] | 1253 | }
|
---|
[2430] | 1254 | gt = g;
|
---|
| 1255 | ga = Math.Abs(gt);
|
---|
[2563] | 1256 | if( (double)(ga)==(double)(0) )
|
---|
[2154] | 1257 | {
|
---|
| 1258 |
|
---|
| 1259 | //
|
---|
[2430] | 1260 | // Diagonal matrix
|
---|
[2154] | 1261 | //
|
---|
[2430] | 1262 | ssmin = ha;
|
---|
| 1263 | ssmax = fa;
|
---|
| 1264 | clt = 1;
|
---|
| 1265 | crt = 1;
|
---|
| 1266 | slt = 0;
|
---|
| 1267 | srt = 0;
|
---|
| 1268 | }
|
---|
| 1269 | else
|
---|
| 1270 | {
|
---|
| 1271 | gasmal = true;
|
---|
[2563] | 1272 | if( (double)(ga)>(double)(fa) )
|
---|
[2154] | 1273 | {
|
---|
[2430] | 1274 | pmax = 2;
|
---|
[2563] | 1275 | if( (double)(fa/ga)<(double)(AP.Math.MachineEpsilon) )
|
---|
[2430] | 1276 | {
|
---|
| 1277 |
|
---|
| 1278 | //
|
---|
| 1279 | // Case of very large GA
|
---|
| 1280 | //
|
---|
| 1281 | gasmal = false;
|
---|
| 1282 | ssmax = ga;
|
---|
[2563] | 1283 | if( (double)(ha)>(double)(1) )
|
---|
[2430] | 1284 | {
|
---|
| 1285 | v = ga/ha;
|
---|
| 1286 | ssmin = fa/v;
|
---|
| 1287 | }
|
---|
| 1288 | else
|
---|
| 1289 | {
|
---|
| 1290 | v = fa/ga;
|
---|
| 1291 | ssmin = v*ha;
|
---|
| 1292 | }
|
---|
| 1293 | clt = 1;
|
---|
| 1294 | slt = ht/gt;
|
---|
| 1295 | srt = 1;
|
---|
| 1296 | crt = ft/gt;
|
---|
| 1297 | }
|
---|
[2154] | 1298 | }
|
---|
[2430] | 1299 | if( gasmal )
|
---|
[2154] | 1300 | {
|
---|
| 1301 |
|
---|
| 1302 | //
|
---|
[2430] | 1303 | // Normal case
|
---|
[2154] | 1304 | //
|
---|
[2430] | 1305 | d = fa-ha;
|
---|
[2563] | 1306 | if( (double)(d)==(double)(fa) )
|
---|
[2430] | 1307 | {
|
---|
| 1308 | l = 1;
|
---|
| 1309 | }
|
---|
| 1310 | else
|
---|
| 1311 | {
|
---|
| 1312 | l = d/fa;
|
---|
| 1313 | }
|
---|
| 1314 | m = gt/ft;
|
---|
| 1315 | t = 2-l;
|
---|
| 1316 | mm = m*m;
|
---|
| 1317 | tt = t*t;
|
---|
| 1318 | s = Math.Sqrt(tt+mm);
|
---|
[2563] | 1319 | if( (double)(l)==(double)(0) )
|
---|
[2154] | 1320 | {
|
---|
[2430] | 1321 | r = Math.Abs(m);
|
---|
[2154] | 1322 | }
|
---|
| 1323 | else
|
---|
| 1324 | {
|
---|
[2430] | 1325 | r = Math.Sqrt(l*l+mm);
|
---|
[2154] | 1326 | }
|
---|
[2430] | 1327 | a = 0.5*(s+r);
|
---|
| 1328 | ssmin = ha/a;
|
---|
| 1329 | ssmax = fa*a;
|
---|
[2563] | 1330 | if( (double)(mm)==(double)(0) )
|
---|
[2430] | 1331 | {
|
---|
| 1332 |
|
---|
| 1333 | //
|
---|
| 1334 | // Note that M is very tiny
|
---|
| 1335 | //
|
---|
[2563] | 1336 | if( (double)(l)==(double)(0) )
|
---|
[2430] | 1337 | {
|
---|
| 1338 | t = extsignbdsqr(2, ft)*extsignbdsqr(1, gt);
|
---|
| 1339 | }
|
---|
| 1340 | else
|
---|
| 1341 | {
|
---|
| 1342 | t = gt/extsignbdsqr(d, ft)+m/t;
|
---|
| 1343 | }
|
---|
| 1344 | }
|
---|
| 1345 | else
|
---|
| 1346 | {
|
---|
| 1347 | t = (m/(s+t)+m/(r+l))*(1+a);
|
---|
| 1348 | }
|
---|
| 1349 | l = Math.Sqrt(t*t+4);
|
---|
| 1350 | crt = 2/l;
|
---|
| 1351 | srt = t/l;
|
---|
| 1352 | clt = (crt+srt*m)/a;
|
---|
| 1353 | v = ht/ft;
|
---|
| 1354 | slt = v*srt/a;
|
---|
[2154] | 1355 | }
|
---|
| 1356 | }
|
---|
[2430] | 1357 | if( swp )
|
---|
| 1358 | {
|
---|
| 1359 | csl = srt;
|
---|
| 1360 | snl = crt;
|
---|
| 1361 | csr = slt;
|
---|
| 1362 | snr = clt;
|
---|
| 1363 | }
|
---|
| 1364 | else
|
---|
| 1365 | {
|
---|
| 1366 | csl = clt;
|
---|
| 1367 | snl = slt;
|
---|
| 1368 | csr = crt;
|
---|
| 1369 | snr = srt;
|
---|
| 1370 | }
|
---|
| 1371 |
|
---|
| 1372 | //
|
---|
| 1373 | // Correct signs of SSMAX and SSMIN
|
---|
| 1374 | //
|
---|
| 1375 | if( pmax==1 )
|
---|
| 1376 | {
|
---|
| 1377 | tsign = extsignbdsqr(1, csr)*extsignbdsqr(1, csl)*extsignbdsqr(1, f);
|
---|
| 1378 | }
|
---|
| 1379 | if( pmax==2 )
|
---|
| 1380 | {
|
---|
| 1381 | tsign = extsignbdsqr(1, snr)*extsignbdsqr(1, csl)*extsignbdsqr(1, g);
|
---|
| 1382 | }
|
---|
| 1383 | if( pmax==3 )
|
---|
| 1384 | {
|
---|
| 1385 | tsign = extsignbdsqr(1, snr)*extsignbdsqr(1, snl)*extsignbdsqr(1, h);
|
---|
| 1386 | }
|
---|
| 1387 | ssmax = extsignbdsqr(ssmax, tsign);
|
---|
| 1388 | ssmin = extsignbdsqr(ssmin, tsign*extsignbdsqr(1, f)*extsignbdsqr(1, h));
|
---|
[2154] | 1389 | }
|
---|
| 1390 | }
|
---|
| 1391 | }
|
---|