Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/estnorm.cs

Last change on this file was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 12.5 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace 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}
Note: See TracBrowser for help on using the repository browser.