[2563] | 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 estnorm
|
---|
| 32 | {
|
---|
| 33 | /*************************************************************************
|
---|
| 34 | Matrix norm estimation
|
---|
| 35 |
|
---|
| 36 | The algorithm estimates the 1-norm of square matrix A on the assumption
|
---|
| 37 | that the multiplication of matrix A by the vector is available (the
|
---|
| 38 | iterative method is used). It is recommended to use this algorithm if it
|
---|
| 39 | is hard to calculate matrix elements explicitly (for example, when
|
---|
| 40 | estimating the inverse matrix norm).
|
---|
| 41 |
|
---|
| 42 | The algorithm uses back communication for multiplying the vector by the
|
---|
| 43 | matrix. If KASE=0 after returning from a subroutine, its execution was
|
---|
| 44 | completed successfully, otherwise it is required to multiply the returned
|
---|
| 45 | vector by matrix A and call the subroutine again.
|
---|
| 46 |
|
---|
| 47 | The DemoIterativeEstimateNorm subroutine shows a simple example.
|
---|
| 48 |
|
---|
| 49 | Parameters:
|
---|
| 50 | N - size of matrix A.
|
---|
| 51 | V - vector. It is initialized by the subroutine on the first
|
---|
| 52 | call. It is then passed into it on repeated calls.
|
---|
| 53 | X - if KASE<>0, it contains the vector to be replaced by:
|
---|
| 54 | A * X, if KASE=1
|
---|
| 55 | A^T * X, if KASE=2
|
---|
| 56 | Array whose index ranges within [1..N].
|
---|
| 57 | ISGN - vector. It is initialized by the subroutine on the first
|
---|
| 58 | call. It is then passed into it on repeated calls.
|
---|
| 59 | EST - if KASE=0, it contains the lower boundary of the matrix
|
---|
| 60 | norm estimate.
|
---|
| 61 | KASE - on the first call, it should be equal to 0. After the last
|
---|
| 62 | return, it is equal to 0 (EST contains the matrix norm),
|
---|
| 63 | on intermediate returns it can be equal to 1 or 2 depending
|
---|
| 64 | on the operation to be performed on vector X.
|
---|
| 65 |
|
---|
| 66 | -- LAPACK auxiliary routine (version 3.0) --
|
---|
| 67 | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
---|
| 68 | Courant Institute, Argonne National Lab, and Rice University
|
---|
| 69 | February 29, 1992
|
---|
| 70 | *************************************************************************/
|
---|
| 71 | public static void iterativeestimate1norm(int n,
|
---|
| 72 | ref double[] v,
|
---|
| 73 | ref double[] x,
|
---|
| 74 | ref int[] isgn,
|
---|
| 75 | ref double est,
|
---|
| 76 | ref int kase)
|
---|
| 77 | {
|
---|
| 78 | int itmax = 0;
|
---|
| 79 | int i = 0;
|
---|
| 80 | double t = 0;
|
---|
| 81 | bool flg = new bool();
|
---|
| 82 | int positer = 0;
|
---|
| 83 | int posj = 0;
|
---|
| 84 | int posjlast = 0;
|
---|
| 85 | int posjump = 0;
|
---|
| 86 | int posaltsgn = 0;
|
---|
| 87 | int posestold = 0;
|
---|
| 88 | int postemp = 0;
|
---|
| 89 | int i_ = 0;
|
---|
| 90 |
|
---|
| 91 | itmax = 5;
|
---|
| 92 | posaltsgn = n+1;
|
---|
| 93 | posestold = n+2;
|
---|
| 94 | postemp = n+3;
|
---|
| 95 | positer = n+1;
|
---|
| 96 | posj = n+2;
|
---|
| 97 | posjlast = n+3;
|
---|
| 98 | posjump = n+4;
|
---|
| 99 | if( kase==0 )
|
---|
| 100 | {
|
---|
| 101 | v = new double[n+3+1];
|
---|
| 102 | x = new double[n+1];
|
---|
| 103 | isgn = new int[n+4+1];
|
---|
| 104 | t = (double)(1)/(double)(n);
|
---|
| 105 | for(i=1; i<=n; i++)
|
---|
| 106 | {
|
---|
| 107 | x[i] = t;
|
---|
| 108 | }
|
---|
| 109 | kase = 1;
|
---|
| 110 | isgn[posjump] = 1;
|
---|
| 111 | return;
|
---|
| 112 | }
|
---|
| 113 |
|
---|
| 114 | //
|
---|
| 115 | // ................ ENTRY (JUMP = 1)
|
---|
| 116 | // FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 117 | //
|
---|
| 118 | if( isgn[posjump]==1 )
|
---|
| 119 | {
|
---|
| 120 | if( n==1 )
|
---|
| 121 | {
|
---|
| 122 | v[1] = x[1];
|
---|
| 123 | est = Math.Abs(v[1]);
|
---|
| 124 | kase = 0;
|
---|
| 125 | return;
|
---|
| 126 | }
|
---|
| 127 | est = 0;
|
---|
| 128 | for(i=1; i<=n; i++)
|
---|
| 129 | {
|
---|
| 130 | est = est+Math.Abs(x[i]);
|
---|
| 131 | }
|
---|
| 132 | for(i=1; i<=n; i++)
|
---|
| 133 | {
|
---|
| 134 | if( (double)(x[i])>=(double)(0) )
|
---|
| 135 | {
|
---|
| 136 | x[i] = 1;
|
---|
| 137 | }
|
---|
| 138 | else
|
---|
| 139 | {
|
---|
| 140 | x[i] = -1;
|
---|
| 141 | }
|
---|
| 142 | isgn[i] = Math.Sign(x[i]);
|
---|
| 143 | }
|
---|
| 144 | kase = 2;
|
---|
| 145 | isgn[posjump] = 2;
|
---|
| 146 | return;
|
---|
| 147 | }
|
---|
| 148 |
|
---|
| 149 | //
|
---|
| 150 | // ................ ENTRY (JUMP = 2)
|
---|
| 151 | // FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
|
---|
| 152 | //
|
---|
| 153 | if( isgn[posjump]==2 )
|
---|
| 154 | {
|
---|
| 155 | isgn[posj] = 1;
|
---|
| 156 | for(i=2; i<=n; i++)
|
---|
| 157 | {
|
---|
| 158 | if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
|
---|
| 159 | {
|
---|
| 160 | isgn[posj] = i;
|
---|
| 161 | }
|
---|
| 162 | }
|
---|
| 163 | isgn[positer] = 2;
|
---|
| 164 |
|
---|
| 165 | //
|
---|
| 166 | // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
|
---|
| 167 | //
|
---|
| 168 | for(i=1; i<=n; i++)
|
---|
| 169 | {
|
---|
| 170 | x[i] = 0;
|
---|
| 171 | }
|
---|
| 172 | x[isgn[posj]] = 1;
|
---|
| 173 | kase = 1;
|
---|
| 174 | isgn[posjump] = 3;
|
---|
| 175 | return;
|
---|
| 176 | }
|
---|
| 177 |
|
---|
| 178 | //
|
---|
| 179 | // ................ ENTRY (JUMP = 3)
|
---|
| 180 | // X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 181 | //
|
---|
| 182 | if( isgn[posjump]==3 )
|
---|
| 183 | {
|
---|
| 184 | for(i_=1; i_<=n;i_++)
|
---|
| 185 | {
|
---|
| 186 | v[i_] = x[i_];
|
---|
| 187 | }
|
---|
| 188 | v[posestold] = est;
|
---|
| 189 | est = 0;
|
---|
| 190 | for(i=1; i<=n; i++)
|
---|
| 191 | {
|
---|
| 192 | est = est+Math.Abs(v[i]);
|
---|
| 193 | }
|
---|
| 194 | flg = false;
|
---|
| 195 | for(i=1; i<=n; i++)
|
---|
| 196 | {
|
---|
| 197 | if( (double)(x[i])>=(double)(0) & isgn[i]<0 | (double)(x[i])<(double)(0) & isgn[i]>=0 )
|
---|
| 198 | {
|
---|
| 199 | flg = true;
|
---|
| 200 | }
|
---|
| 201 | }
|
---|
| 202 |
|
---|
| 203 | //
|
---|
| 204 | // REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
|
---|
| 205 | // OR MAY BE CYCLING.
|
---|
| 206 | //
|
---|
| 207 | if( !flg | (double)(est)<=(double)(v[posestold]) )
|
---|
| 208 | {
|
---|
| 209 | v[posaltsgn] = 1;
|
---|
| 210 | for(i=1; i<=n; i++)
|
---|
| 211 | {
|
---|
| 212 | x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
|
---|
| 213 | v[posaltsgn] = -v[posaltsgn];
|
---|
| 214 | }
|
---|
| 215 | kase = 1;
|
---|
| 216 | isgn[posjump] = 5;
|
---|
| 217 | return;
|
---|
| 218 | }
|
---|
| 219 | for(i=1; i<=n; i++)
|
---|
| 220 | {
|
---|
| 221 | if( (double)(x[i])>=(double)(0) )
|
---|
| 222 | {
|
---|
| 223 | x[i] = 1;
|
---|
| 224 | isgn[i] = 1;
|
---|
| 225 | }
|
---|
| 226 | else
|
---|
| 227 | {
|
---|
| 228 | x[i] = -1;
|
---|
| 229 | isgn[i] = -1;
|
---|
| 230 | }
|
---|
| 231 | }
|
---|
| 232 | kase = 2;
|
---|
| 233 | isgn[posjump] = 4;
|
---|
| 234 | return;
|
---|
| 235 | }
|
---|
| 236 |
|
---|
| 237 | //
|
---|
| 238 | // ................ ENTRY (JUMP = 4)
|
---|
| 239 | // X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
|
---|
| 240 | //
|
---|
| 241 | if( isgn[posjump]==4 )
|
---|
| 242 | {
|
---|
| 243 | isgn[posjlast] = isgn[posj];
|
---|
| 244 | isgn[posj] = 1;
|
---|
| 245 | for(i=2; i<=n; i++)
|
---|
| 246 | {
|
---|
| 247 | if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
|
---|
| 248 | {
|
---|
| 249 | isgn[posj] = i;
|
---|
| 250 | }
|
---|
| 251 | }
|
---|
| 252 | if( (double)(x[isgn[posjlast]])!=(double)(Math.Abs(x[isgn[posj]])) & isgn[positer]<itmax )
|
---|
| 253 | {
|
---|
| 254 | isgn[positer] = isgn[positer]+1;
|
---|
| 255 | for(i=1; i<=n; i++)
|
---|
| 256 | {
|
---|
| 257 | x[i] = 0;
|
---|
| 258 | }
|
---|
| 259 | x[isgn[posj]] = 1;
|
---|
| 260 | kase = 1;
|
---|
| 261 | isgn[posjump] = 3;
|
---|
| 262 | return;
|
---|
| 263 | }
|
---|
| 264 |
|
---|
| 265 | //
|
---|
| 266 | // ITERATION COMPLETE. FINAL STAGE.
|
---|
| 267 | //
|
---|
| 268 | v[posaltsgn] = 1;
|
---|
| 269 | for(i=1; i<=n; i++)
|
---|
| 270 | {
|
---|
| 271 | x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
|
---|
| 272 | v[posaltsgn] = -v[posaltsgn];
|
---|
| 273 | }
|
---|
| 274 | kase = 1;
|
---|
| 275 | isgn[posjump] = 5;
|
---|
| 276 | return;
|
---|
| 277 | }
|
---|
| 278 |
|
---|
| 279 | //
|
---|
| 280 | // ................ ENTRY (JUMP = 5)
|
---|
| 281 | // X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 282 | //
|
---|
| 283 | if( isgn[posjump]==5 )
|
---|
| 284 | {
|
---|
| 285 | v[postemp] = 0;
|
---|
| 286 | for(i=1; i<=n; i++)
|
---|
| 287 | {
|
---|
| 288 | v[postemp] = v[postemp]+Math.Abs(x[i]);
|
---|
| 289 | }
|
---|
| 290 | v[postemp] = 2*v[postemp]/(3*n);
|
---|
| 291 | if( (double)(v[postemp])>(double)(est) )
|
---|
| 292 | {
|
---|
| 293 | for(i_=1; i_<=n;i_++)
|
---|
| 294 | {
|
---|
| 295 | v[i_] = x[i_];
|
---|
| 296 | }
|
---|
| 297 | est = v[postemp];
|
---|
| 298 | }
|
---|
| 299 | kase = 0;
|
---|
| 300 | return;
|
---|
| 301 | }
|
---|
| 302 | }
|
---|
| 303 |
|
---|
| 304 |
|
---|
| 305 | /*************************************************************************
|
---|
| 306 | Example of usage of an IterativeEstimateNorm subroutine
|
---|
| 307 |
|
---|
| 308 | Input parameters:
|
---|
| 309 | A - matrix.
|
---|
| 310 | Array whose indexes range within [1..N, 1..N].
|
---|
| 311 |
|
---|
| 312 | Return:
|
---|
| 313 | Matrix norm estimated by the subroutine.
|
---|
| 314 |
|
---|
| 315 | -- ALGLIB --
|
---|
| 316 | Copyright 2005 by Bochkanov Sergey
|
---|
| 317 | *************************************************************************/
|
---|
| 318 | public static double demoiterativeestimate1norm(ref double[,] a,
|
---|
| 319 | int n)
|
---|
| 320 | {
|
---|
| 321 | double result = 0;
|
---|
| 322 | int i = 0;
|
---|
| 323 | double s = 0;
|
---|
| 324 | double[] x = new double[0];
|
---|
| 325 | double[] t = new double[0];
|
---|
| 326 | double[] v = new double[0];
|
---|
| 327 | int[] iv = new int[0];
|
---|
| 328 | int kase = 0;
|
---|
| 329 | int i_ = 0;
|
---|
| 330 |
|
---|
| 331 | kase = 0;
|
---|
| 332 | t = new double[n+1];
|
---|
| 333 | iterativeestimate1norm(n, ref v, ref x, ref iv, ref result, ref kase);
|
---|
| 334 | while( kase!=0 )
|
---|
| 335 | {
|
---|
| 336 | if( kase==1 )
|
---|
| 337 | {
|
---|
| 338 | for(i=1; i<=n; i++)
|
---|
| 339 | {
|
---|
| 340 | s = 0.0;
|
---|
| 341 | for(i_=1; i_<=n;i_++)
|
---|
| 342 | {
|
---|
| 343 | s += a[i,i_]*x[i_];
|
---|
| 344 | }
|
---|
| 345 | t[i] = s;
|
---|
| 346 | }
|
---|
| 347 | }
|
---|
| 348 | else
|
---|
| 349 | {
|
---|
| 350 | for(i=1; i<=n; i++)
|
---|
| 351 | {
|
---|
| 352 | s = 0.0;
|
---|
| 353 | for(i_=1; i_<=n;i_++)
|
---|
| 354 | {
|
---|
| 355 | s += a[i_,i]*x[i_];
|
---|
| 356 | }
|
---|
| 357 | t[i] = s;
|
---|
| 358 | }
|
---|
| 359 | }
|
---|
| 360 | for(i_=1; i_<=n;i_++)
|
---|
| 361 | {
|
---|
| 362 | x[i_] = t[i_];
|
---|
| 363 | }
|
---|
| 364 | iterativeestimate1norm(n, ref v, ref x, ref iv, ref result, ref kase);
|
---|
| 365 | }
|
---|
| 366 | return result;
|
---|
| 367 | }
|
---|
| 368 | }
|
---|
| 369 | }
|
---|