Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/tdbisinv.cs @ 2645

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

extracted external libraries and adapted dependent plugins (ticket #837)
File size: 90.7 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 tdbisinv
32    {
33        /*************************************************************************
34        Subroutine for finding the tridiagonal matrix eigenvalues/vectors in a
35        given half-interval (A, B] by using bisection and inverse iteration.
36
37        Input parameters:
38            D       -   the main diagonal of a tridiagonal matrix.
39                        Array whose index ranges within [0..N-1].
40            E       -   the secondary diagonal of a tridiagonal matrix.
41                        Array whose index ranges within [0..N-2].
42            N       -   size of matrix, N>=0.
43            ZNeeded -   flag controlling whether the eigenvectors are needed or not.
44                        If ZNeeded is equal to:
45                         * 0, the eigenvectors are not needed;
46                         * 1, the eigenvectors of a tridiagonal matrix are multiplied
47                           by the square matrix Z. It is used if the tridiagonal
48                           matrix is obtained by the similarity transformation
49                           of a symmetric matrix.
50                         * 2, the eigenvectors of a tridiagonal matrix replace matrix Z.
51            A, B    -   half-interval (A, B] to search eigenvalues in.
52            Z       -   if ZNeeded is equal to:
53                         * 0, Z isn't used and remains unchanged;
54                         * 1, Z contains the square matrix (array whose indexes range
55                           within [0..N-1, 0..N-1]) which reduces the given symmetric
56                           matrix to tridiagonal form;
57                         * 2, Z isn't used (but changed on the exit).
58
59        Output parameters:
60            D       -   array of the eigenvalues found.
61                        Array whose index ranges within [0..M-1].
62            M       -   number of eigenvalues found in the given half-interval (M>=0).
63            Z       -   if ZNeeded is equal to:
64                         * 0, doesn't contain any information;
65                         * 1, contains the product of a given NxN matrix Z (from the
66                           left) and NxM matrix of the eigenvectors found (from the
67                           right). Array whose indexes range within [0..N-1, 0..M-1].
68                         * 2, contains the matrix of the eigenvectors found.
69                           Array whose indexes range within [0..N-1, 0..M-1].
70
71        Result:
72
73            True, if successful. In that case, M contains the number of eigenvalues
74            in the given half-interval (could be equal to 0), D contains the eigenvalues,
75            Z contains the eigenvectors (if needed).
76            It should be noted that the subroutine changes the size of arrays D and Z.
77
78            False, if the bisection method subroutine wasn't able to find the
79            eigenvalues in the given interval or if the inverse iteration subroutine
80            wasn't able to find all the corresponding eigenvectors. In that case,
81            the eigenvalues and eigenvectors are not returned, M is equal to 0.
82
83          -- ALGLIB --
84             Copyright 31.03.2008 by Bochkanov Sergey
85        *************************************************************************/
86        public static bool smatrixtdevdr(ref double[] d,
87            ref double[] e,
88            int n,
89            int zneeded,
90            double a,
91            double b,
92            ref int m,
93            ref double[,] z)
94        {
95            bool result = new bool();
96            int errorcode = 0;
97            int nsplit = 0;
98            int i = 0;
99            int j = 0;
100            int k = 0;
101            int cr = 0;
102            int[] iblock = new int[0];
103            int[] isplit = new int[0];
104            int[] ifail = new int[0];
105            double[] d1 = new double[0];
106            double[] e1 = new double[0];
107            double[] w = new double[0];
108            double[,] z2 = new double[0,0];
109            double[,] z3 = new double[0,0];
110            double v = 0;
111            int i_ = 0;
112            int i1_ = 0;
113
114            System.Diagnostics.Debug.Assert(zneeded>=0 & zneeded<=2, "SMatrixTDEVDR: incorrect ZNeeded!");
115           
116            //
117            // Special cases
118            //
119            if( (double)(b)<=(double)(a) )
120            {
121                m = 0;
122                result = true;
123                return result;
124            }
125            if( n<=0 )
126            {
127                m = 0;
128                result = true;
129                return result;
130            }
131           
132            //
133            // Copy D,E to D1, E1
134            //
135            d1 = new double[n+1];
136            i1_ = (0) - (1);
137            for(i_=1; i_<=n;i_++)
138            {
139                d1[i_] = d[i_+i1_];
140            }
141            if( n>1 )
142            {
143                e1 = new double[n-1+1];
144                i1_ = (0) - (1);
145                for(i_=1; i_<=n-1;i_++)
146                {
147                    e1[i_] = e[i_+i1_];
148                }
149            }
150           
151            //
152            // No eigen vectors
153            //
154            if( zneeded==0 )
155            {
156                result = internalbisectioneigenvalues(d1, e1, n, 2, 1, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
157                if( !result | m==0 )
158                {
159                    m = 0;
160                    return result;
161                }
162                d = new double[m-1+1];
163                i1_ = (1) - (0);
164                for(i_=0; i_<=m-1;i_++)
165                {
166                    d[i_] = w[i_+i1_];
167                }
168                return result;
169            }
170           
171            //
172            // Eigen vectors are multiplied by Z
173            //
174            if( zneeded==1 )
175            {
176               
177                //
178                // Find eigen pairs
179                //
180                result = internalbisectioneigenvalues(d1, e1, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
181                if( !result | m==0 )
182                {
183                    m = 0;
184                    return result;
185                }
186                internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
187                if( cr!=0 )
188                {
189                    m = 0;
190                    result = false;
191                    return result;
192                }
193               
194                //
195                // Sort eigen values and vectors
196                //
197                for(i=1; i<=m; i++)
198                {
199                    k = i;
200                    for(j=i; j<=m; j++)
201                    {
202                        if( (double)(w[j])<(double)(w[k]) )
203                        {
204                            k = j;
205                        }
206                    }
207                    v = w[i];
208                    w[i] = w[k];
209                    w[k] = v;
210                    for(j=1; j<=n; j++)
211                    {
212                        v = z2[j,i];
213                        z2[j,i] = z2[j,k];
214                        z2[j,k] = v;
215                    }
216                }
217               
218                //
219                // Transform Z2 and overwrite Z
220                //
221                z3 = new double[m+1, n+1];
222                for(i=1; i<=m; i++)
223                {
224                    for(i_=1; i_<=n;i_++)
225                    {
226                        z3[i,i_] = z2[i_,i];
227                    }
228                }
229                for(i=1; i<=n; i++)
230                {
231                    for(j=1; j<=m; j++)
232                    {
233                        i1_ = (1)-(0);
234                        v = 0.0;
235                        for(i_=0; i_<=n-1;i_++)
236                        {
237                            v += z[i-1,i_]*z3[j,i_+i1_];
238                        }
239                        z2[i,j] = v;
240                    }
241                }
242                z = new double[n-1+1, m-1+1];
243                for(i=1; i<=m; i++)
244                {
245                    i1_ = (1) - (0);
246                    for(i_=0; i_<=n-1;i_++)
247                    {
248                        z[i_,i-1] = z2[i_+i1_,i];
249                    }
250                }
251               
252                //
253                // Store W
254                //
255                d = new double[m-1+1];
256                for(i=1; i<=m; i++)
257                {
258                    d[i-1] = w[i];
259                }
260                return result;
261            }
262           
263            //
264            // Eigen vectors are stored in Z
265            //
266            if( zneeded==2 )
267            {
268               
269                //
270                // Find eigen pairs
271                //
272                result = internalbisectioneigenvalues(d1, e1, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
273                if( !result | m==0 )
274                {
275                    m = 0;
276                    return result;
277                }
278                internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
279                if( cr!=0 )
280                {
281                    m = 0;
282                    result = false;
283                    return result;
284                }
285               
286                //
287                // Sort eigen values and vectors
288                //
289                for(i=1; i<=m; i++)
290                {
291                    k = i;
292                    for(j=i; j<=m; j++)
293                    {
294                        if( (double)(w[j])<(double)(w[k]) )
295                        {
296                            k = j;
297                        }
298                    }
299                    v = w[i];
300                    w[i] = w[k];
301                    w[k] = v;
302                    for(j=1; j<=n; j++)
303                    {
304                        v = z2[j,i];
305                        z2[j,i] = z2[j,k];
306                        z2[j,k] = v;
307                    }
308                }
309               
310                //
311                // Store W
312                //
313                d = new double[m-1+1];
314                for(i=1; i<=m; i++)
315                {
316                    d[i-1] = w[i];
317                }
318                z = new double[n-1+1, m-1+1];
319                for(i=1; i<=m; i++)
320                {
321                    i1_ = (1) - (0);
322                    for(i_=0; i_<=n-1;i_++)
323                    {
324                        z[i_,i-1] = z2[i_+i1_,i];
325                    }
326                }
327                return result;
328            }
329            result = false;
330            return result;
331        }
332
333
334        /*************************************************************************
335        Subroutine for finding tridiagonal matrix eigenvalues/vectors with given
336        indexes (in ascending order) by using the bisection and inverse iteraion.
337
338        Input parameters:
339            D       -   the main diagonal of a tridiagonal matrix.
340                        Array whose index ranges within [0..N-1].
341            E       -   the secondary diagonal of a tridiagonal matrix.
342                        Array whose index ranges within [0..N-2].
343            N       -   size of matrix. N>=0.
344            ZNeeded -   flag controlling whether the eigenvectors are needed or not.
345                        If ZNeeded is equal to:
346                         * 0, the eigenvectors are not needed;
347                         * 1, the eigenvectors of a tridiagonal matrix are multiplied
348                           by the square matrix Z. It is used if the
349                           tridiagonal matrix is obtained by the similarity transformation
350                           of a symmetric matrix.
351                         * 2, the eigenvectors of a tridiagonal matrix replace
352                           matrix Z.
353            I1, I2  -   index interval for searching (from I1 to I2).
354                        0 <= I1 <= I2 <= N-1.
355            Z       -   if ZNeeded is equal to:
356                         * 0, Z isn't used and remains unchanged;
357                         * 1, Z contains the square matrix (array whose indexes range within [0..N-1, 0..N-1])
358                           which reduces the given symmetric matrix to  tridiagonal form;
359                         * 2, Z isn't used (but changed on the exit).
360
361        Output parameters:
362            D       -   array of the eigenvalues found.
363                        Array whose index ranges within [0..I2-I1].
364            Z       -   if ZNeeded is equal to:
365                         * 0, doesn't contain any information;
366                         * 1, contains the product of a given NxN matrix Z (from the left) and
367                           Nx(I2-I1) matrix of the eigenvectors found (from the right).
368                           Array whose indexes range within [0..N-1, 0..I2-I1].
369                         * 2, contains the matrix of the eigenvalues found.
370                           Array whose indexes range within [0..N-1, 0..I2-I1].
371
372
373        Result:
374
375            True, if successful. In that case, D contains the eigenvalues,
376            Z contains the eigenvectors (if needed).
377            It should be noted that the subroutine changes the size of arrays D and Z.
378
379            False, if the bisection method subroutine wasn't able to find the eigenvalues
380            in the given interval or if the inverse iteration subroutine wasn't able
381            to find all the corresponding eigenvectors. In that case, the eigenvalues
382            and eigenvectors are not returned.
383
384          -- ALGLIB --
385             Copyright 25.12.2005 by Bochkanov Sergey
386        *************************************************************************/
387        public static bool smatrixtdevdi(ref double[] d,
388            ref double[] e,
389            int n,
390            int zneeded,
391            int i1,
392            int i2,
393            ref double[,] z)
394        {
395            bool result = new bool();
396            int errorcode = 0;
397            int nsplit = 0;
398            int i = 0;
399            int j = 0;
400            int k = 0;
401            int m = 0;
402            int cr = 0;
403            int[] iblock = new int[0];
404            int[] isplit = new int[0];
405            int[] ifail = new int[0];
406            double[] w = new double[0];
407            double[] d1 = new double[0];
408            double[] e1 = new double[0];
409            double[,] z2 = new double[0,0];
410            double[,] z3 = new double[0,0];
411            double v = 0;
412            int i_ = 0;
413            int i1_ = 0;
414
415            System.Diagnostics.Debug.Assert(0<=i1 & i1<=i2 & i2<n, "SMatrixTDEVDI: incorrect I1/I2!");
416           
417            //
418            // Copy D,E to D1, E1
419            //
420            d1 = new double[n+1];
421            i1_ = (0) - (1);
422            for(i_=1; i_<=n;i_++)
423            {
424                d1[i_] = d[i_+i1_];
425            }
426            if( n>1 )
427            {
428                e1 = new double[n-1+1];
429                i1_ = (0) - (1);
430                for(i_=1; i_<=n-1;i_++)
431                {
432                    e1[i_] = e[i_+i1_];
433                }
434            }
435           
436            //
437            // No eigen vectors
438            //
439            if( zneeded==0 )
440            {
441                result = internalbisectioneigenvalues(d1, e1, n, 3, 1, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
442                if( !result )
443                {
444                    return result;
445                }
446                if( m!=i2-i1+1 )
447                {
448                    result = false;
449                    return result;
450                }
451                d = new double[m-1+1];
452                for(i=1; i<=m; i++)
453                {
454                    d[i-1] = w[i];
455                }
456                return result;
457            }
458           
459            //
460            // Eigen vectors are multiplied by Z
461            //
462            if( zneeded==1 )
463            {
464               
465                //
466                // Find eigen pairs
467                //
468                result = internalbisectioneigenvalues(d1, e1, n, 3, 2, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
469                if( !result )
470                {
471                    return result;
472                }
473                if( m!=i2-i1+1 )
474                {
475                    result = false;
476                    return result;
477                }
478                internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
479                if( cr!=0 )
480                {
481                    result = false;
482                    return result;
483                }
484               
485                //
486                // Sort eigen values and vectors
487                //
488                for(i=1; i<=m; i++)
489                {
490                    k = i;
491                    for(j=i; j<=m; j++)
492                    {
493                        if( (double)(w[j])<(double)(w[k]) )
494                        {
495                            k = j;
496                        }
497                    }
498                    v = w[i];
499                    w[i] = w[k];
500                    w[k] = v;
501                    for(j=1; j<=n; j++)
502                    {
503                        v = z2[j,i];
504                        z2[j,i] = z2[j,k];
505                        z2[j,k] = v;
506                    }
507                }
508               
509                //
510                // Transform Z2 and overwrite Z
511                //
512                z3 = new double[m+1, n+1];
513                for(i=1; i<=m; i++)
514                {
515                    for(i_=1; i_<=n;i_++)
516                    {
517                        z3[i,i_] = z2[i_,i];
518                    }
519                }
520                for(i=1; i<=n; i++)
521                {
522                    for(j=1; j<=m; j++)
523                    {
524                        i1_ = (1)-(0);
525                        v = 0.0;
526                        for(i_=0; i_<=n-1;i_++)
527                        {
528                            v += z[i-1,i_]*z3[j,i_+i1_];
529                        }
530                        z2[i,j] = v;
531                    }
532                }
533                z = new double[n-1+1, m-1+1];
534                for(i=1; i<=m; i++)
535                {
536                    i1_ = (1) - (0);
537                    for(i_=0; i_<=n-1;i_++)
538                    {
539                        z[i_,i-1] = z2[i_+i1_,i];
540                    }
541                }
542               
543                //
544                // Store W
545                //
546                d = new double[m-1+1];
547                for(i=1; i<=m; i++)
548                {
549                    d[i-1] = w[i];
550                }
551                return result;
552            }
553           
554            //
555            // Eigen vectors are stored in Z
556            //
557            if( zneeded==2 )
558            {
559               
560                //
561                // Find eigen pairs
562                //
563                result = internalbisectioneigenvalues(d1, e1, n, 3, 2, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
564                if( !result )
565                {
566                    return result;
567                }
568                if( m!=i2-i1+1 )
569                {
570                    result = false;
571                    return result;
572                }
573                internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
574                if( cr!=0 )
575                {
576                    result = false;
577                    return result;
578                }
579               
580                //
581                // Sort eigen values and vectors
582                //
583                for(i=1; i<=m; i++)
584                {
585                    k = i;
586                    for(j=i; j<=m; j++)
587                    {
588                        if( (double)(w[j])<(double)(w[k]) )
589                        {
590                            k = j;
591                        }
592                    }
593                    v = w[i];
594                    w[i] = w[k];
595                    w[k] = v;
596                    for(j=1; j<=n; j++)
597                    {
598                        v = z2[j,i];
599                        z2[j,i] = z2[j,k];
600                        z2[j,k] = v;
601                    }
602                }
603               
604                //
605                // Store Z
606                //
607                z = new double[n-1+1, m-1+1];
608                for(i=1; i<=m; i++)
609                {
610                    i1_ = (1) - (0);
611                    for(i_=0; i_<=n-1;i_++)
612                    {
613                        z[i_,i-1] = z2[i_+i1_,i];
614                    }
615                }
616               
617                //
618                // Store W
619                //
620                d = new double[m-1+1];
621                for(i=1; i<=m; i++)
622                {
623                    d[i-1] = w[i];
624                }
625                return result;
626            }
627            result = false;
628            return result;
629        }
630
631
632        public static bool tridiagonaleigenvaluesandvectorsininterval(ref double[] d,
633            ref double[] e,
634            int n,
635            int zneeded,
636            double a,
637            double b,
638            ref int m,
639            ref double[,] z)
640        {
641            bool result = new bool();
642            int errorcode = 0;
643            int nsplit = 0;
644            int i = 0;
645            int j = 0;
646            int k = 0;
647            int cr = 0;
648            int[] iblock = new int[0];
649            int[] isplit = new int[0];
650            int[] ifail = new int[0];
651            double[] w = new double[0];
652            double[,] z2 = new double[0,0];
653            double[,] z3 = new double[0,0];
654            double v = 0;
655            int i_ = 0;
656
657           
658            //
659            // No eigen vectors
660            //
661            if( zneeded==0 )
662            {
663                result = internalbisectioneigenvalues(d, e, n, 2, 1, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
664                if( !result | m==0 )
665                {
666                    m = 0;
667                    return result;
668                }
669                d = new double[m+1];
670                for(i=1; i<=m; i++)
671                {
672                    d[i] = w[i];
673                }
674                return result;
675            }
676           
677            //
678            // Eigen vectors are multiplied by Z
679            //
680            if( zneeded==1 )
681            {
682               
683                //
684                // Find eigen pairs
685                //
686                result = internalbisectioneigenvalues(d, e, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
687                if( !result | m==0 )
688                {
689                    m = 0;
690                    return result;
691                }
692                internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
693                if( cr!=0 )
694                {
695                    m = 0;
696                    result = false;
697                    return result;
698                }
699               
700                //
701                // Sort eigen values and vectors
702                //
703                for(i=1; i<=m; i++)
704                {
705                    k = i;
706                    for(j=i; j<=m; j++)
707                    {
708                        if( (double)(w[j])<(double)(w[k]) )
709                        {
710                            k = j;
711                        }
712                    }
713                    v = w[i];
714                    w[i] = w[k];
715                    w[k] = v;
716                    for(j=1; j<=n; j++)
717                    {
718                        v = z2[j,i];
719                        z2[j,i] = z2[j,k];
720                        z2[j,k] = v;
721                    }
722                }
723               
724                //
725                // Transform Z2 and overwrite Z
726                //
727                z3 = new double[m+1, n+1];
728                for(i=1; i<=m; i++)
729                {
730                    for(i_=1; i_<=n;i_++)
731                    {
732                        z3[i,i_] = z2[i_,i];
733                    }
734                }
735                for(i=1; i<=n; i++)
736                {
737                    for(j=1; j<=m; j++)
738                    {
739                        v = 0.0;
740                        for(i_=1; i_<=n;i_++)
741                        {
742                            v += z[i,i_]*z3[j,i_];
743                        }
744                        z2[i,j] = v;
745                    }
746                }
747                z = new double[n+1, m+1];
748                for(i=1; i<=m; i++)
749                {
750                    for(i_=1; i_<=n;i_++)
751                    {
752                        z[i_,i] = z2[i_,i];
753                    }
754                }
755               
756                //
757                // Store W
758                //
759                d = new double[m+1];
760                for(i=1; i<=m; i++)
761                {
762                    d[i] = w[i];
763                }
764                return result;
765            }
766           
767            //
768            // Eigen vectors are stored in Z
769            //
770            if( zneeded==2 )
771            {
772               
773                //
774                // Find eigen pairs
775                //
776                result = internalbisectioneigenvalues(d, e, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
777                if( !result | m==0 )
778                {
779                    m = 0;
780                    return result;
781                }
782                internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z, ref ifail, ref cr);
783                if( cr!=0 )
784                {
785                    m = 0;
786                    result = false;
787                    return result;
788                }
789               
790                //
791                // Sort eigen values and vectors
792                //
793                for(i=1; i<=m; i++)
794                {
795                    k = i;
796                    for(j=i; j<=m; j++)
797                    {
798                        if( (double)(w[j])<(double)(w[k]) )
799                        {
800                            k = j;
801                        }
802                    }
803                    v = w[i];
804                    w[i] = w[k];
805                    w[k] = v;
806                    for(j=1; j<=n; j++)
807                    {
808                        v = z[j,i];
809                        z[j,i] = z[j,k];
810                        z[j,k] = v;
811                    }
812                }
813               
814                //
815                // Store W
816                //
817                d = new double[m+1];
818                for(i=1; i<=m; i++)
819                {
820                    d[i] = w[i];
821                }
822                return result;
823            }
824            result = false;
825            return result;
826        }
827
828
829        public static bool tridiagonaleigenvaluesandvectorsbyindexes(ref double[] d,
830            ref double[] e,
831            int n,
832            int zneeded,
833            int i1,
834            int i2,
835            ref double[,] z)
836        {
837            bool result = new bool();
838            int errorcode = 0;
839            int nsplit = 0;
840            int i = 0;
841            int j = 0;
842            int k = 0;
843            int m = 0;
844            int cr = 0;
845            int[] iblock = new int[0];
846            int[] isplit = new int[0];
847            int[] ifail = new int[0];
848            double[] w = new double[0];
849            double[,] z2 = new double[0,0];
850            double[,] z3 = new double[0,0];
851            double v = 0;
852            int i_ = 0;
853
854           
855            //
856            // No eigen vectors
857            //
858            if( zneeded==0 )
859            {
860                result = internalbisectioneigenvalues(d, e, n, 3, 1, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
861                if( !result )
862                {
863                    return result;
864                }
865                d = new double[m+1];
866                for(i=1; i<=m; i++)
867                {
868                    d[i] = w[i];
869                }
870                return result;
871            }
872           
873            //
874            // Eigen vectors are multiplied by Z
875            //
876            if( zneeded==1 )
877            {
878               
879                //
880                // Find eigen pairs
881                //
882                result = internalbisectioneigenvalues(d, e, n, 3, 2, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
883                if( !result )
884                {
885                    return result;
886                }
887                internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
888                if( cr!=0 )
889                {
890                    result = false;
891                    return result;
892                }
893               
894                //
895                // Sort eigen values and vectors
896                //
897                for(i=1; i<=m; i++)
898                {
899                    k = i;
900                    for(j=i; j<=m; j++)
901                    {
902                        if( (double)(w[j])<(double)(w[k]) )
903                        {
904                            k = j;
905                        }
906                    }
907                    v = w[i];
908                    w[i] = w[k];
909                    w[k] = v;
910                    for(j=1; j<=n; j++)
911                    {
912                        v = z2[j,i];
913                        z2[j,i] = z2[j,k];
914                        z2[j,k] = v;
915                    }
916                }
917               
918                //
919                // Transform Z2 and overwrite Z
920                //
921                z3 = new double[m+1, n+1];
922                for(i=1; i<=m; i++)
923                {
924                    for(i_=1; i_<=n;i_++)
925                    {
926                        z3[i,i_] = z2[i_,i];
927                    }
928                }
929                for(i=1; i<=n; i++)
930                {
931                    for(j=1; j<=m; j++)
932                    {
933                        v = 0.0;
934                        for(i_=1; i_<=n;i_++)
935                        {
936                            v += z[i,i_]*z3[j,i_];
937                        }
938                        z2[i,j] = v;
939                    }
940                }
941                z = new double[n+1, m+1];
942                for(i=1; i<=m; i++)
943                {
944                    for(i_=1; i_<=n;i_++)
945                    {
946                        z[i_,i] = z2[i_,i];
947                    }
948                }
949               
950                //
951                // Store W
952                //
953                d = new double[m+1];
954                for(i=1; i<=m; i++)
955                {
956                    d[i] = w[i];
957                }
958                return result;
959            }
960           
961            //
962            // Eigen vectors are stored in Z
963            //
964            if( zneeded==2 )
965            {
966               
967                //
968                // Find eigen pairs
969                //
970                result = internalbisectioneigenvalues(d, e, n, 3, 2, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
971                if( !result )
972                {
973                    return result;
974                }
975                internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z, ref ifail, ref cr);
976                if( cr!=0 )
977                {
978                    result = false;
979                    return result;
980                }
981               
982                //
983                // Sort eigen values and vectors
984                //
985                for(i=1; i<=m; i++)
986                {
987                    k = i;
988                    for(j=i; j<=m; j++)
989                    {
990                        if( (double)(w[j])<(double)(w[k]) )
991                        {
992                            k = j;
993                        }
994                    }
995                    v = w[i];
996                    w[i] = w[k];
997                    w[k] = v;
998                    for(j=1; j<=n; j++)
999                    {
1000                        v = z[j,i];
1001                        z[j,i] = z[j,k];
1002                        z[j,k] = v;
1003                    }
1004                }
1005               
1006                //
1007                // Store W
1008                //
1009                d = new double[m+1];
1010                for(i=1; i<=m; i++)
1011                {
1012                    d[i] = w[i];
1013                }
1014                return result;
1015            }
1016            result = false;
1017            return result;
1018        }
1019
1020
1021        public static bool internalbisectioneigenvalues(double[] d,
1022            double[] e,
1023            int n,
1024            int irange,
1025            int iorder,
1026            double vl,
1027            double vu,
1028            int il,
1029            int iu,
1030            double abstol,
1031            ref double[] w,
1032            ref int m,
1033            ref int nsplit,
1034            ref int[] iblock,
1035            ref int[] isplit,
1036            ref int errorcode)
1037        {
1038            bool result = new bool();
1039            double fudge = 0;
1040            double relfac = 0;
1041            bool ncnvrg = new bool();
1042            bool toofew = new bool();
1043            int ib = 0;
1044            int ibegin = 0;
1045            int idiscl = 0;
1046            int idiscu = 0;
1047            int ie = 0;
1048            int iend = 0;
1049            int iinfo = 0;
1050            int im = 0;
1051            int iin = 0;
1052            int ioff = 0;
1053            int iout = 0;
1054            int itmax = 0;
1055            int iw = 0;
1056            int iwoff = 0;
1057            int j = 0;
1058            int itmp1 = 0;
1059            int jb = 0;
1060            int jdisc = 0;
1061            int je = 0;
1062            int nwl = 0;
1063            int nwu = 0;
1064            double atoli = 0;
1065            double bnorm = 0;
1066            double gl = 0;
1067            double gu = 0;
1068            double pivmin = 0;
1069            double rtoli = 0;
1070            double safemn = 0;
1071            double tmp1 = 0;
1072            double tmp2 = 0;
1073            double tnorm = 0;
1074            double ulp = 0;
1075            double wkill = 0;
1076            double wl = 0;
1077            double wlu = 0;
1078            double wu = 0;
1079            double wul = 0;
1080            double scalefactor = 0;
1081            double t = 0;
1082            int[] idumma = new int[0];
1083            double[] work = new double[0];
1084            int[] iwork = new int[0];
1085            int[] ia1s2 = new int[0];
1086            double[] ra1s2 = new double[0];
1087            double[,] ra1s2x2 = new double[0,0];
1088            int[,] ia1s2x2 = new int[0,0];
1089            double[] ra1siin = new double[0];
1090            double[] ra2siin = new double[0];
1091            double[] ra3siin = new double[0];
1092            double[] ra4siin = new double[0];
1093            double[,] ra1siinx2 = new double[0,0];
1094            int[,] ia1siinx2 = new int[0,0];
1095            int[] iworkspace = new int[0];
1096            double[] rworkspace = new double[0];
1097            int tmpi = 0;
1098
1099            d = (double[])d.Clone();
1100            e = (double[])e.Clone();
1101
1102           
1103            //
1104            // Quick return if possible
1105            //
1106            m = 0;
1107            if( n==0 )
1108            {
1109                result = true;
1110                return result;
1111            }
1112           
1113            //
1114            // Get machine constants
1115            // NB is the minimum vector length for vector bisection, or 0
1116            // if only scalar is to be done.
1117            //
1118            fudge = 2;
1119            relfac = 2;
1120            safemn = AP.Math.MinRealNumber;
1121            ulp = 2*AP.Math.MachineEpsilon;
1122            rtoli = ulp*relfac;
1123            idumma = new int[1+1];
1124            work = new double[4*n+1];
1125            iwork = new int[3*n+1];
1126            w = new double[n+1];
1127            iblock = new int[n+1];
1128            isplit = new int[n+1];
1129            ia1s2 = new int[2+1];
1130            ra1s2 = new double[2+1];
1131            ra1s2x2 = new double[2+1, 2+1];
1132            ia1s2x2 = new int[2+1, 2+1];
1133            ra1siin = new double[n+1];
1134            ra2siin = new double[n+1];
1135            ra3siin = new double[n+1];
1136            ra4siin = new double[n+1];
1137            ra1siinx2 = new double[n+1, 2+1];
1138            ia1siinx2 = new int[n+1, 2+1];
1139            iworkspace = new int[n+1];
1140            rworkspace = new double[n+1];
1141           
1142            //
1143            // Check for Errors
1144            //
1145            result = false;
1146            errorcode = 0;
1147            if( irange<=0 | irange>=4 )
1148            {
1149                errorcode = -4;
1150            }
1151            if( iorder<=0 | iorder>=3 )
1152            {
1153                errorcode = -5;
1154            }
1155            if( n<0 )
1156            {
1157                errorcode = -3;
1158            }
1159            if( irange==2 & (double)(vl)>=(double)(vu) )
1160            {
1161                errorcode = -6;
1162            }
1163            if( irange==3 & (il<1 | il>Math.Max(1, n)) )
1164            {
1165                errorcode = -8;
1166            }
1167            if( irange==3 & (iu<Math.Min(n, il) | iu>n) )
1168            {
1169                errorcode = -9;
1170            }
1171            if( errorcode!=0 )
1172            {
1173                return result;
1174            }
1175           
1176            //
1177            // Initialize error flags
1178            //
1179            ncnvrg = false;
1180            toofew = false;
1181           
1182            //
1183            // Simplifications:
1184            //
1185            if( irange==3 & il==1 & iu==n )
1186            {
1187                irange = 1;
1188            }
1189           
1190            //
1191            // Special Case when N=1
1192            //
1193            if( n==1 )
1194            {
1195                nsplit = 1;
1196                isplit[1] = 1;
1197                if( irange==2 & ((double)(vl)>=(double)(d[1]) | (double)(vu)<(double)(d[1])) )
1198                {
1199                    m = 0;
1200                }
1201                else
1202                {
1203                    w[1] = d[1];
1204                    iblock[1] = 1;
1205                    m = 1;
1206                }
1207                result = true;
1208                return result;
1209            }
1210           
1211            //
1212            // Scaling
1213            //
1214            t = Math.Abs(d[n]);
1215            for(j=1; j<=n-1; j++)
1216            {
1217                t = Math.Max(t, Math.Abs(d[j]));
1218                t = Math.Max(t, Math.Abs(e[j]));
1219            }
1220            scalefactor = 1;
1221            if( (double)(t)!=(double)(0) )
1222            {
1223                if( (double)(t)>(double)(Math.Sqrt(Math.Sqrt(AP.Math.MinRealNumber))*Math.Sqrt(AP.Math.MaxRealNumber)) )
1224                {
1225                    scalefactor = t;
1226                }
1227                if( (double)(t)<(double)(Math.Sqrt(Math.Sqrt(AP.Math.MaxRealNumber))*Math.Sqrt(AP.Math.MinRealNumber)) )
1228                {
1229                    scalefactor = t;
1230                }
1231                for(j=1; j<=n-1; j++)
1232                {
1233                    d[j] = d[j]/scalefactor;
1234                    e[j] = e[j]/scalefactor;
1235                }
1236                d[n] = d[n]/scalefactor;
1237            }
1238           
1239            //
1240            // Compute Splitting Points
1241            //
1242            nsplit = 1;
1243            work[n] = 0;
1244            pivmin = 1;
1245            for(j=2; j<=n; j++)
1246            {
1247                tmp1 = AP.Math.Sqr(e[j-1]);
1248                if( (double)(Math.Abs(d[j]*d[j-1])*AP.Math.Sqr(ulp)+safemn)>(double)(tmp1) )
1249                {
1250                    isplit[nsplit] = j-1;
1251                    nsplit = nsplit+1;
1252                    work[j-1] = 0;
1253                }
1254                else
1255                {
1256                    work[j-1] = tmp1;
1257                    pivmin = Math.Max(pivmin, tmp1);
1258                }
1259            }
1260            isplit[nsplit] = n;
1261            pivmin = pivmin*safemn;
1262           
1263            //
1264            // Compute Interval and ATOLI
1265            //
1266            if( irange==3 )
1267            {
1268               
1269                //
1270                // RANGE='I': Compute the interval containing eigenvalues
1271                //     IL through IU.
1272                //
1273                // Compute Gershgorin interval for entire (split) matrix
1274                // and use it as the initial interval
1275                //
1276                gu = d[1];
1277                gl = d[1];
1278                tmp1 = 0;
1279                for(j=1; j<=n-1; j++)
1280                {
1281                    tmp2 = Math.Sqrt(work[j]);
1282                    gu = Math.Max(gu, d[j]+tmp1+tmp2);
1283                    gl = Math.Min(gl, d[j]-tmp1-tmp2);
1284                    tmp1 = tmp2;
1285                }
1286                gu = Math.Max(gu, d[n]+tmp1);
1287                gl = Math.Min(gl, d[n]-tmp1);
1288                tnorm = Math.Max(Math.Abs(gl), Math.Abs(gu));
1289                gl = gl-fudge*tnorm*ulp*n-fudge*2*pivmin;
1290                gu = gu+fudge*tnorm*ulp*n+fudge*pivmin;
1291               
1292                //
1293                // Compute Iteration parameters
1294                //
1295                itmax = (int)Math.Ceiling((Math.Log(tnorm+pivmin)-Math.Log(pivmin))/Math.Log(2))+2;
1296                if( (double)(abstol)<=(double)(0) )
1297                {
1298                    atoli = ulp*tnorm;
1299                }
1300                else
1301                {
1302                    atoli = abstol;
1303                }
1304                work[n+1] = gl;
1305                work[n+2] = gl;
1306                work[n+3] = gu;
1307                work[n+4] = gu;
1308                work[n+5] = gl;
1309                work[n+6] = gu;
1310                iwork[1] = -1;
1311                iwork[2] = -1;
1312                iwork[3] = n+1;
1313                iwork[4] = n+1;
1314                iwork[5] = il-1;
1315                iwork[6] = iu;
1316               
1317                //
1318                // Calling DLAEBZ
1319                //
1320                // DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
1321                //    WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
1322                //    IWORK, W, IBLOCK, IINFO )
1323                //
1324                ia1s2[1] = iwork[5];
1325                ia1s2[2] = iwork[6];
1326                ra1s2[1] = work[n+5];
1327                ra1s2[2] = work[n+6];
1328                ra1s2x2[1,1] = work[n+1];
1329                ra1s2x2[2,1] = work[n+2];
1330                ra1s2x2[1,2] = work[n+3];
1331                ra1s2x2[2,2] = work[n+4];
1332                ia1s2x2[1,1] = iwork[1];
1333                ia1s2x2[2,1] = iwork[2];
1334                ia1s2x2[1,2] = iwork[3];
1335                ia1s2x2[2,2] = iwork[4];
1336                internaldlaebz(3, itmax, n, 2, 2, atoli, rtoli, pivmin, ref d, ref e, ref work, ref ia1s2, ref ra1s2x2, ref ra1s2, ref iout, ref ia1s2x2, ref w, ref iblock, ref iinfo);
1337                iwork[5] = ia1s2[1];
1338                iwork[6] = ia1s2[2];
1339                work[n+5] = ra1s2[1];
1340                work[n+6] = ra1s2[2];
1341                work[n+1] = ra1s2x2[1,1];
1342                work[n+2] = ra1s2x2[2,1];
1343                work[n+3] = ra1s2x2[1,2];
1344                work[n+4] = ra1s2x2[2,2];
1345                iwork[1] = ia1s2x2[1,1];
1346                iwork[2] = ia1s2x2[2,1];
1347                iwork[3] = ia1s2x2[1,2];
1348                iwork[4] = ia1s2x2[2,2];
1349                if( iwork[6]==iu )
1350                {
1351                    wl = work[n+1];
1352                    wlu = work[n+3];
1353                    nwl = iwork[1];
1354                    wu = work[n+4];
1355                    wul = work[n+2];
1356                    nwu = iwork[4];
1357                }
1358                else
1359                {
1360                    wl = work[n+2];
1361                    wlu = work[n+4];
1362                    nwl = iwork[2];
1363                    wu = work[n+3];
1364                    wul = work[n+1];
1365                    nwu = iwork[3];
1366                }
1367                if( nwl<0 | nwl>=n | nwu<1 | nwu>n )
1368                {
1369                    errorcode = 4;
1370                    result = false;
1371                    return result;
1372                }
1373            }
1374            else
1375            {
1376               
1377                //
1378                // RANGE='A' or 'V' -- Set ATOLI
1379                //
1380                tnorm = Math.Max(Math.Abs(d[1])+Math.Abs(e[1]), Math.Abs(d[n])+Math.Abs(e[n-1]));
1381                for(j=2; j<=n-1; j++)
1382                {
1383                    tnorm = Math.Max(tnorm, Math.Abs(d[j])+Math.Abs(e[j-1])+Math.Abs(e[j]));
1384                }
1385                if( (double)(abstol)<=(double)(0) )
1386                {
1387                    atoli = ulp*tnorm;
1388                }
1389                else
1390                {
1391                    atoli = abstol;
1392                }
1393                if( irange==2 )
1394                {
1395                    wl = vl;
1396                    wu = vu;
1397                }
1398                else
1399                {
1400                    wl = 0;
1401                    wu = 0;
1402                }
1403            }
1404           
1405            //
1406            // Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
1407            // NWL accumulates the number of eigenvalues .le. WL,
1408            // NWU accumulates the number of eigenvalues .le. WU
1409            //
1410            m = 0;
1411            iend = 0;
1412            errorcode = 0;
1413            nwl = 0;
1414            nwu = 0;
1415            for(jb=1; jb<=nsplit; jb++)
1416            {
1417                ioff = iend;
1418                ibegin = ioff+1;
1419                iend = isplit[jb];
1420                iin = iend-ioff;
1421                if( iin==1 )
1422                {
1423                   
1424                    //
1425                    // Special Case -- IIN=1
1426                    //
1427                    if( irange==1 | (double)(wl)>=(double)(d[ibegin]-pivmin) )
1428                    {
1429                        nwl = nwl+1;
1430                    }
1431                    if( irange==1 | (double)(wu)>=(double)(d[ibegin]-pivmin) )
1432                    {
1433                        nwu = nwu+1;
1434                    }
1435                    if( irange==1 | (double)(wl)<(double)(d[ibegin]-pivmin) & (double)(wu)>=(double)(d[ibegin]-pivmin) )
1436                    {
1437                        m = m+1;
1438                        w[m] = d[ibegin];
1439                        iblock[m] = jb;
1440                    }
1441                }
1442                else
1443                {
1444                   
1445                    //
1446                    // General Case -- IIN > 1
1447                    //
1448                    // Compute Gershgorin Interval
1449                    // and use it as the initial interval
1450                    //
1451                    gu = d[ibegin];
1452                    gl = d[ibegin];
1453                    tmp1 = 0;
1454                    for(j=ibegin; j<=iend-1; j++)
1455                    {
1456                        tmp2 = Math.Abs(e[j]);
1457                        gu = Math.Max(gu, d[j]+tmp1+tmp2);
1458                        gl = Math.Min(gl, d[j]-tmp1-tmp2);
1459                        tmp1 = tmp2;
1460                    }
1461                    gu = Math.Max(gu, d[iend]+tmp1);
1462                    gl = Math.Min(gl, d[iend]-tmp1);
1463                    bnorm = Math.Max(Math.Abs(gl), Math.Abs(gu));
1464                    gl = gl-fudge*bnorm*ulp*iin-fudge*pivmin;
1465                    gu = gu+fudge*bnorm*ulp*iin+fudge*pivmin;
1466                   
1467                    //
1468                    // Compute ATOLI for the current submatrix
1469                    //
1470                    if( (double)(abstol)<=(double)(0) )
1471                    {
1472                        atoli = ulp*Math.Max(Math.Abs(gl), Math.Abs(gu));
1473                    }
1474                    else
1475                    {
1476                        atoli = abstol;
1477                    }
1478                    if( irange>1 )
1479                    {
1480                        if( (double)(gu)<(double)(wl) )
1481                        {
1482                            nwl = nwl+iin;
1483                            nwu = nwu+iin;
1484                            continue;
1485                        }
1486                        gl = Math.Max(gl, wl);
1487                        gu = Math.Min(gu, wu);
1488                        if( (double)(gl)>=(double)(gu) )
1489                        {
1490                            continue;
1491                        }
1492                    }
1493                   
1494                    //
1495                    // Set Up Initial Interval
1496                    //
1497                    work[n+1] = gl;
1498                    work[n+iin+1] = gu;
1499                   
1500                    //
1501                    // Calling DLAEBZ
1502                    //
1503                    // CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
1504                    //    D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
1505                    //    IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
1506                    //    IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
1507                    //
1508                    for(tmpi=1; tmpi<=iin; tmpi++)
1509                    {
1510                        ra1siin[tmpi] = d[ibegin-1+tmpi];
1511                        if( ibegin-1+tmpi<n )
1512                        {
1513                            ra2siin[tmpi] = e[ibegin-1+tmpi];
1514                        }
1515                        ra3siin[tmpi] = work[ibegin-1+tmpi];
1516                        ra1siinx2[tmpi,1] = work[n+tmpi];
1517                        ra1siinx2[tmpi,2] = work[n+tmpi+iin];
1518                        ra4siin[tmpi] = work[n+2*iin+tmpi];
1519                        rworkspace[tmpi] = w[m+tmpi];
1520                        iworkspace[tmpi] = iblock[m+tmpi];
1521                        ia1siinx2[tmpi,1] = iwork[tmpi];
1522                        ia1siinx2[tmpi,2] = iwork[tmpi+iin];
1523                    }
1524                    internaldlaebz(1, 0, iin, iin, 1, atoli, rtoli, pivmin, ref ra1siin, ref ra2siin, ref ra3siin, ref idumma, ref ra1siinx2, ref ra4siin, ref im, ref ia1siinx2, ref rworkspace, ref iworkspace, ref iinfo);
1525                    for(tmpi=1; tmpi<=iin; tmpi++)
1526                    {
1527                        work[n+tmpi] = ra1siinx2[tmpi,1];
1528                        work[n+tmpi+iin] = ra1siinx2[tmpi,2];
1529                        work[n+2*iin+tmpi] = ra4siin[tmpi];
1530                        w[m+tmpi] = rworkspace[tmpi];
1531                        iblock[m+tmpi] = iworkspace[tmpi];
1532                        iwork[tmpi] = ia1siinx2[tmpi,1];
1533                        iwork[tmpi+iin] = ia1siinx2[tmpi,2];
1534                    }
1535                    nwl = nwl+iwork[1];
1536                    nwu = nwu+iwork[iin+1];
1537                    iwoff = m-iwork[1];
1538                   
1539                    //
1540                    // Compute Eigenvalues
1541                    //
1542                    itmax = (int)Math.Ceiling((Math.Log(gu-gl+pivmin)-Math.Log(pivmin))/Math.Log(2))+2;
1543                   
1544                    //
1545                    // Calling DLAEBZ
1546                    //
1547                    //CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
1548                    //    D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
1549                    //    IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
1550                    //    IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
1551                    //
1552                    for(tmpi=1; tmpi<=iin; tmpi++)
1553                    {
1554                        ra1siin[tmpi] = d[ibegin-1+tmpi];
1555                        if( ibegin-1+tmpi<n )
1556                        {
1557                            ra2siin[tmpi] = e[ibegin-1+tmpi];
1558                        }
1559                        ra3siin[tmpi] = work[ibegin-1+tmpi];
1560                        ra1siinx2[tmpi,1] = work[n+tmpi];
1561                        ra1siinx2[tmpi,2] = work[n+tmpi+iin];
1562                        ra4siin[tmpi] = work[n+2*iin+tmpi];
1563                        rworkspace[tmpi] = w[m+tmpi];
1564                        iworkspace[tmpi] = iblock[m+tmpi];
1565                        ia1siinx2[tmpi,1] = iwork[tmpi];
1566                        ia1siinx2[tmpi,2] = iwork[tmpi+iin];
1567                    }
1568                    internaldlaebz(2, itmax, iin, iin, 1, atoli, rtoli, pivmin, ref ra1siin, ref ra2siin, ref ra3siin, ref idumma, ref ra1siinx2, ref ra4siin, ref iout, ref ia1siinx2, ref rworkspace, ref iworkspace, ref iinfo);
1569                    for(tmpi=1; tmpi<=iin; tmpi++)
1570                    {
1571                        work[n+tmpi] = ra1siinx2[tmpi,1];
1572                        work[n+tmpi+iin] = ra1siinx2[tmpi,2];
1573                        work[n+2*iin+tmpi] = ra4siin[tmpi];
1574                        w[m+tmpi] = rworkspace[tmpi];
1575                        iblock[m+tmpi] = iworkspace[tmpi];
1576                        iwork[tmpi] = ia1siinx2[tmpi,1];
1577                        iwork[tmpi+iin] = ia1siinx2[tmpi,2];
1578                    }
1579                   
1580                    //
1581                    // Copy Eigenvalues Into W and IBLOCK
1582                    // Use -JB for block number for unconverged eigenvalues.
1583                    //
1584                    for(j=1; j<=iout; j++)
1585                    {
1586                        tmp1 = 0.5*(work[j+n]+work[j+iin+n]);
1587                       
1588                        //
1589                        // Flag non-convergence.
1590                        //
1591                        if( j>iout-iinfo )
1592                        {
1593                            ncnvrg = true;
1594                            ib = -jb;
1595                        }
1596                        else
1597                        {
1598                            ib = jb;
1599                        }
1600                        for(je=iwork[j]+1+iwoff; je<=iwork[j+iin]+iwoff; je++)
1601                        {
1602                            w[je] = tmp1;
1603                            iblock[je] = ib;
1604                        }
1605                    }
1606                    m = m+im;
1607                }
1608            }
1609           
1610            //
1611            // If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
1612            // If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
1613            //
1614            if( irange==3 )
1615            {
1616                im = 0;
1617                idiscl = il-1-nwl;
1618                idiscu = nwu-iu;
1619                if( idiscl>0 | idiscu>0 )
1620                {
1621                    for(je=1; je<=m; je++)
1622                    {
1623                        if( (double)(w[je])<=(double)(wlu) & idiscl>0 )
1624                        {
1625                            idiscl = idiscl-1;
1626                        }
1627                        else
1628                        {
1629                            if( (double)(w[je])>=(double)(wul) & idiscu>0 )
1630                            {
1631                                idiscu = idiscu-1;
1632                            }
1633                            else
1634                            {
1635                                im = im+1;
1636                                w[im] = w[je];
1637                                iblock[im] = iblock[je];
1638                            }
1639                        }
1640                    }
1641                    m = im;
1642                }
1643                if( idiscl>0 | idiscu>0 )
1644                {
1645                   
1646                    //
1647                    // Code to deal with effects of bad arithmetic:
1648                    // Some low eigenvalues to be discarded are not in (WL,WLU],
1649                    // or high eigenvalues to be discarded are not in (WUL,WU]
1650                    // so just kill off the smallest IDISCL/largest IDISCU
1651                    // eigenvalues, by simply finding the smallest/largest
1652                    // eigenvalue(s).
1653                    //
1654                    // (If N(w) is monotone non-decreasing, this should never
1655                    //  happen.)
1656                    //
1657                    if( idiscl>0 )
1658                    {
1659                        wkill = wu;
1660                        for(jdisc=1; jdisc<=idiscl; jdisc++)
1661                        {
1662                            iw = 0;
1663                            for(je=1; je<=m; je++)
1664                            {
1665                                if( iblock[je]!=0 & ((double)(w[je])<(double)(wkill) | iw==0) )
1666                                {
1667                                    iw = je;
1668                                    wkill = w[je];
1669                                }
1670                            }
1671                            iblock[iw] = 0;
1672                        }
1673                    }
1674                    if( idiscu>0 )
1675                    {
1676                        wkill = wl;
1677                        for(jdisc=1; jdisc<=idiscu; jdisc++)
1678                        {
1679                            iw = 0;
1680                            for(je=1; je<=m; je++)
1681                            {
1682                                if( iblock[je]!=0 & ((double)(w[je])>(double)(wkill) | iw==0) )
1683                                {
1684                                    iw = je;
1685                                    wkill = w[je];
1686                                }
1687                            }
1688                            iblock[iw] = 0;
1689                        }
1690                    }
1691                    im = 0;
1692                    for(je=1; je<=m; je++)
1693                    {
1694                        if( iblock[je]!=0 )
1695                        {
1696                            im = im+1;
1697                            w[im] = w[je];
1698                            iblock[im] = iblock[je];
1699                        }
1700                    }
1701                    m = im;
1702                }
1703                if( idiscl<0 | idiscu<0 )
1704                {
1705                    toofew = true;
1706                }
1707            }
1708           
1709            //
1710            // If ORDER='B', do nothing -- the eigenvalues are already sorted
1711            //    by block.
1712            // If ORDER='E', sort the eigenvalues from smallest to largest
1713            //
1714            if( iorder==1 & nsplit>1 )
1715            {
1716                for(je=1; je<=m-1; je++)
1717                {
1718                    ie = 0;
1719                    tmp1 = w[je];
1720                    for(j=je+1; j<=m; j++)
1721                    {
1722                        if( (double)(w[j])<(double)(tmp1) )
1723                        {
1724                            ie = j;
1725                            tmp1 = w[j];
1726                        }
1727                    }
1728                    if( ie!=0 )
1729                    {
1730                        itmp1 = iblock[ie];
1731                        w[ie] = w[je];
1732                        iblock[ie] = iblock[je];
1733                        w[je] = tmp1;
1734                        iblock[je] = itmp1;
1735                    }
1736                }
1737            }
1738            for(j=1; j<=m; j++)
1739            {
1740                w[j] = w[j]*scalefactor;
1741            }
1742            errorcode = 0;
1743            if( ncnvrg )
1744            {
1745                errorcode = errorcode+1;
1746            }
1747            if( toofew )
1748            {
1749                errorcode = errorcode+2;
1750            }
1751            result = errorcode==0;
1752            return result;
1753        }
1754
1755
1756        public static void internaldstein(int n,
1757            ref double[] d,
1758            double[] e,
1759            int m,
1760            double[] w,
1761            ref int[] iblock,
1762            ref int[] isplit,
1763            ref double[,] z,
1764            ref int[] ifail,
1765            ref int info)
1766        {
1767            int maxits = 0;
1768            int extra = 0;
1769            int b1 = 0;
1770            int blksiz = 0;
1771            int bn = 0;
1772            int gpind = 0;
1773            int i = 0;
1774            int iinfo = 0;
1775            int its = 0;
1776            int j = 0;
1777            int j1 = 0;
1778            int jblk = 0;
1779            int jmax = 0;
1780            int nblk = 0;
1781            int nrmchk = 0;
1782            double dtpcrt = 0;
1783            double eps = 0;
1784            double eps1 = 0;
1785            double nrm = 0;
1786            double onenrm = 0;
1787            double ortol = 0;
1788            double pertol = 0;
1789            double scl = 0;
1790            double sep = 0;
1791            double tol = 0;
1792            double xj = 0;
1793            double xjm = 0;
1794            double ztr = 0;
1795            double[] work1 = new double[0];
1796            double[] work2 = new double[0];
1797            double[] work3 = new double[0];
1798            double[] work4 = new double[0];
1799            double[] work5 = new double[0];
1800            int[] iwork = new int[0];
1801            bool tmpcriterion = new bool();
1802            int ti = 0;
1803            int i1 = 0;
1804            int i2 = 0;
1805            double v = 0;
1806            int i_ = 0;
1807            int i1_ = 0;
1808
1809            e = (double[])e.Clone();
1810            w = (double[])w.Clone();
1811
1812            maxits = 5;
1813            extra = 2;
1814            work1 = new double[Math.Max(n, 1)+1];
1815            work2 = new double[Math.Max(n-1, 1)+1];
1816            work3 = new double[Math.Max(n, 1)+1];
1817            work4 = new double[Math.Max(n, 1)+1];
1818            work5 = new double[Math.Max(n, 1)+1];
1819            iwork = new int[Math.Max(n, 1)+1];
1820            ifail = new int[Math.Max(m, 1)+1];
1821            z = new double[Math.Max(n, 1)+1, Math.Max(m, 1)+1];
1822           
1823            //
1824            // Test the input parameters.
1825            //
1826            info = 0;
1827            for(i=1; i<=m; i++)
1828            {
1829                ifail[i] = 0;
1830            }
1831            if( n<0 )
1832            {
1833                info = -1;
1834                return;
1835            }
1836            if( m<0 | m>n )
1837            {
1838                info = -4;
1839                return;
1840            }
1841            for(j=2; j<=m; j++)
1842            {
1843                if( iblock[j]<iblock[j-1] )
1844                {
1845                    info = -6;
1846                    break;
1847                }
1848                if( iblock[j]==iblock[j-1] & (double)(w[j])<(double)(w[j-1]) )
1849                {
1850                    info = -5;
1851                    break;
1852                }
1853            }
1854            if( info!=0 )
1855            {
1856                return;
1857            }
1858           
1859            //
1860            // Quick return if possible
1861            //
1862            if( n==0 | m==0 )
1863            {
1864                return;
1865            }
1866            if( n==1 )
1867            {
1868                z[1,1] = 1;
1869                return;
1870            }
1871           
1872            //
1873            // Some preparations
1874            //
1875            ti = n-1;
1876            for(i_=1; i_<=ti;i_++)
1877            {
1878                work1[i_] = e[i_];
1879            }
1880            e = new double[n+1];
1881            for(i_=1; i_<=ti;i_++)
1882            {
1883                e[i_] = work1[i_];
1884            }
1885            for(i_=1; i_<=m;i_++)
1886            {
1887                work1[i_] = w[i_];
1888            }
1889            w = new double[n+1];
1890            for(i_=1; i_<=m;i_++)
1891            {
1892                w[i_] = work1[i_];
1893            }
1894           
1895            //
1896            // Get machine constants.
1897            //
1898            eps = AP.Math.MachineEpsilon;
1899           
1900            //
1901            // Compute eigenvectors of matrix blocks.
1902            //
1903            j1 = 1;
1904            for(nblk=1; nblk<=iblock[m]; nblk++)
1905            {
1906               
1907                //
1908                // Find starting and ending indices of block nblk.
1909                //
1910                if( nblk==1 )
1911                {
1912                    b1 = 1;
1913                }
1914                else
1915                {
1916                    b1 = isplit[nblk-1]+1;
1917                }
1918                bn = isplit[nblk];
1919                blksiz = bn-b1+1;
1920                if( blksiz!=1 )
1921                {
1922                   
1923                    //
1924                    // Compute reorthogonalization criterion and stopping criterion.
1925                    //
1926                    gpind = b1;
1927                    onenrm = Math.Abs(d[b1])+Math.Abs(e[b1]);
1928                    onenrm = Math.Max(onenrm, Math.Abs(d[bn])+Math.Abs(e[bn-1]));
1929                    for(i=b1+1; i<=bn-1; i++)
1930                    {
1931                        onenrm = Math.Max(onenrm, Math.Abs(d[i])+Math.Abs(e[i-1])+Math.Abs(e[i]));
1932                    }
1933                    ortol = 0.001*onenrm;
1934                    dtpcrt = Math.Sqrt(0.1/blksiz);
1935                }
1936               
1937                //
1938                // Loop through eigenvalues of block nblk.
1939                //
1940                jblk = 0;
1941                for(j=j1; j<=m; j++)
1942                {
1943                    if( iblock[j]!=nblk )
1944                    {
1945                        j1 = j;
1946                        break;
1947                    }
1948                    jblk = jblk+1;
1949                    xj = w[j];
1950                    if( blksiz==1 )
1951                    {
1952                       
1953                        //
1954                        // Skip all the work if the block size is one.
1955                        //
1956                        work1[1] = 1;
1957                    }
1958                    else
1959                    {
1960                       
1961                        //
1962                        // If eigenvalues j and j-1 are too close, add a relatively
1963                        // small perturbation.
1964                        //
1965                        if( jblk>1 )
1966                        {
1967                            eps1 = Math.Abs(eps*xj);
1968                            pertol = 10*eps1;
1969                            sep = xj-xjm;
1970                            if( (double)(sep)<(double)(pertol) )
1971                            {
1972                                xj = xjm+pertol;
1973                            }
1974                        }
1975                        its = 0;
1976                        nrmchk = 0;
1977                       
1978                        //
1979                        // Get random starting vector.
1980                        //
1981                        for(ti=1; ti<=blksiz; ti++)
1982                        {
1983                            work1[ti] = 2*AP.Math.RandomReal()-1;
1984                        }
1985                       
1986                        //
1987                        // Copy the matrix T so it won't be destroyed in factorization.
1988                        //
1989                        for(ti=1; ti<=blksiz-1; ti++)
1990                        {
1991                            work2[ti] = e[b1+ti-1];
1992                            work3[ti] = e[b1+ti-1];
1993                            work4[ti] = d[b1+ti-1];
1994                        }
1995                        work4[blksiz] = d[b1+blksiz-1];
1996                       
1997                        //
1998                        // Compute LU factors with partial pivoting  ( PT = LU )
1999                        //
2000                        tol = 0;
2001                        tdininternaldlagtf(blksiz, ref work4, xj, ref work2, ref work3, tol, ref work5, ref iwork, ref iinfo);
2002                       
2003                        //
2004                        // Update iteration count.
2005                        //
2006                        do
2007                        {
2008                            its = its+1;
2009                            if( its>maxits )
2010                            {
2011                               
2012                                //
2013                                // If stopping criterion was not satisfied, update info and
2014                                // store eigenvector number in array ifail.
2015                                //
2016                                info = info+1;
2017                                ifail[info] = j;
2018                                break;
2019                            }
2020                           
2021                            //
2022                            // Normalize and scale the righthand side vector Pb.
2023                            //
2024                            v = 0;
2025                            for(ti=1; ti<=blksiz; ti++)
2026                            {
2027                                v = v+Math.Abs(work1[ti]);
2028                            }
2029                            scl = blksiz*onenrm*Math.Max(eps, Math.Abs(work4[blksiz]))/v;
2030                            for(i_=1; i_<=blksiz;i_++)
2031                            {
2032                                work1[i_] = scl*work1[i_];
2033                            }
2034                           
2035                            //
2036                            // Solve the system LU = Pb.
2037                            //
2038                            tdininternaldlagts(blksiz, ref work4, ref work2, ref work3, ref work5, ref iwork, ref work1, ref tol, ref iinfo);
2039                           
2040                            //
2041                            // Reorthogonalize by modified Gram-Schmidt if eigenvalues are
2042                            // close enough.
2043                            //
2044                            if( jblk!=1 )
2045                            {
2046                                if( (double)(Math.Abs(xj-xjm))>(double)(ortol) )
2047                                {
2048                                    gpind = j;
2049                                }
2050                                if( gpind!=j )
2051                                {
2052                                    for(i=gpind; i<=j-1; i++)
2053                                    {
2054                                        i1 = b1;
2055                                        i2 = b1+blksiz-1;
2056                                        i1_ = (i1)-(1);
2057                                        ztr = 0.0;
2058                                        for(i_=1; i_<=blksiz;i_++)
2059                                        {
2060                                            ztr += work1[i_]*z[i_+i1_,i];
2061                                        }
2062                                        i1_ = (i1) - (1);
2063                                        for(i_=1; i_<=blksiz;i_++)
2064                                        {
2065                                            work1[i_] = work1[i_] - ztr*z[i_+i1_,i];
2066                                        }
2067                                    }
2068                                }
2069                            }
2070                           
2071                            //
2072                            // Check the infinity norm of the iterate.
2073                            //
2074                            jmax = blas.vectoridxabsmax(ref work1, 1, blksiz);
2075                            nrm = Math.Abs(work1[jmax]);
2076                           
2077                            //
2078                            // Continue for additional iterations after norm reaches
2079                            // stopping criterion.
2080                            //
2081                            tmpcriterion = false;
2082                            if( (double)(nrm)<(double)(dtpcrt) )
2083                            {
2084                                tmpcriterion = true;
2085                            }
2086                            else
2087                            {
2088                                nrmchk = nrmchk+1;
2089                                if( nrmchk<extra+1 )
2090                                {
2091                                    tmpcriterion = true;
2092                                }
2093                            }
2094                        }
2095                        while( tmpcriterion );
2096                       
2097                        //
2098                        // Accept iterate as jth eigenvector.
2099                        //
2100                        scl = 1/blas.vectornorm2(ref work1, 1, blksiz);
2101                        jmax = blas.vectoridxabsmax(ref work1, 1, blksiz);
2102                        if( (double)(work1[jmax])<(double)(0) )
2103                        {
2104                            scl = -scl;
2105                        }
2106                        for(i_=1; i_<=blksiz;i_++)
2107                        {
2108                            work1[i_] = scl*work1[i_];
2109                        }
2110                    }
2111                    for(i=1; i<=n; i++)
2112                    {
2113                        z[i,j] = 0;
2114                    }
2115                    for(i=1; i<=blksiz; i++)
2116                    {
2117                        z[b1+i-1,j] = work1[i];
2118                    }
2119                   
2120                    //
2121                    // Save the shift to check eigenvalue spacing at next
2122                    // iteration.
2123                    //
2124                    xjm = xj;
2125                }
2126            }
2127        }
2128
2129
2130        private static void tdininternaldlagtf(int n,
2131            ref double[] a,
2132            double lambda,
2133            ref double[] b,
2134            ref double[] c,
2135            double tol,
2136            ref double[] d,
2137            ref int[] iin,
2138            ref int info)
2139        {
2140            int k = 0;
2141            double eps = 0;
2142            double mult = 0;
2143            double piv1 = 0;
2144            double piv2 = 0;
2145            double scale1 = 0;
2146            double scale2 = 0;
2147            double temp = 0;
2148            double tl = 0;
2149
2150            info = 0;
2151            if( n<0 )
2152            {
2153                info = -1;
2154                return;
2155            }
2156            if( n==0 )
2157            {
2158                return;
2159            }
2160            a[1] = a[1]-lambda;
2161            iin[n] = 0;
2162            if( n==1 )
2163            {
2164                if( (double)(a[1])==(double)(0) )
2165                {
2166                    iin[1] = 1;
2167                }
2168                return;
2169            }
2170            eps = AP.Math.MachineEpsilon;
2171            tl = Math.Max(tol, eps);
2172            scale1 = Math.Abs(a[1])+Math.Abs(b[1]);
2173            for(k=1; k<=n-1; k++)
2174            {
2175                a[k+1] = a[k+1]-lambda;
2176                scale2 = Math.Abs(c[k])+Math.Abs(a[k+1]);
2177                if( k<n-1 )
2178                {
2179                    scale2 = scale2+Math.Abs(b[k+1]);
2180                }
2181                if( (double)(a[k])==(double)(0) )
2182                {
2183                    piv1 = 0;
2184                }
2185                else
2186                {
2187                    piv1 = Math.Abs(a[k])/scale1;
2188                }
2189                if( (double)(c[k])==(double)(0) )
2190                {
2191                    iin[k] = 0;
2192                    piv2 = 0;
2193                    scale1 = scale2;
2194                    if( k<n-1 )
2195                    {
2196                        d[k] = 0;
2197                    }
2198                }
2199                else
2200                {
2201                    piv2 = Math.Abs(c[k])/scale2;
2202                    if( (double)(piv2)<=(double)(piv1) )
2203                    {
2204                        iin[k] = 0;
2205                        scale1 = scale2;
2206                        c[k] = c[k]/a[k];
2207                        a[k+1] = a[k+1]-c[k]*b[k];
2208                        if( k<n-1 )
2209                        {
2210                            d[k] = 0;
2211                        }
2212                    }
2213                    else
2214                    {
2215                        iin[k] = 1;
2216                        mult = a[k]/c[k];
2217                        a[k] = c[k];
2218                        temp = a[k+1];
2219                        a[k+1] = b[k]-mult*temp;
2220                        if( k<n-1 )
2221                        {
2222                            d[k] = b[k+1];
2223                            b[k+1] = -(mult*d[k]);
2224                        }
2225                        b[k] = temp;
2226                        c[k] = mult;
2227                    }
2228                }
2229                if( (double)(Math.Max(piv1, piv2))<=(double)(tl) & iin[n]==0 )
2230                {
2231                    iin[n] = k;
2232                }
2233            }
2234            if( (double)(Math.Abs(a[n]))<=(double)(scale1*tl) & iin[n]==0 )
2235            {
2236                iin[n] = n;
2237            }
2238        }
2239
2240
2241        private static void tdininternaldlagts(int n,
2242            ref double[] a,
2243            ref double[] b,
2244            ref double[] c,
2245            ref double[] d,
2246            ref int[] iin,
2247            ref double[] y,
2248            ref double tol,
2249            ref int info)
2250        {
2251            int k = 0;
2252            double absak = 0;
2253            double ak = 0;
2254            double bignum = 0;
2255            double eps = 0;
2256            double pert = 0;
2257            double sfmin = 0;
2258            double temp = 0;
2259
2260            info = 0;
2261            if( n<0 )
2262            {
2263                info = -1;
2264                return;
2265            }
2266            if( n==0 )
2267            {
2268                return;
2269            }
2270            eps = AP.Math.MachineEpsilon;
2271            sfmin = AP.Math.MinRealNumber;
2272            bignum = 1/sfmin;
2273            if( (double)(tol)<=(double)(0) )
2274            {
2275                tol = Math.Abs(a[1]);
2276                if( n>1 )
2277                {
2278                    tol = Math.Max(tol, Math.Max(Math.Abs(a[2]), Math.Abs(b[1])));
2279                }
2280                for(k=3; k<=n; k++)
2281                {
2282                    tol = Math.Max(tol, Math.Max(Math.Abs(a[k]), Math.Max(Math.Abs(b[k-1]), Math.Abs(d[k-2]))));
2283                }
2284                tol = tol*eps;
2285                if( (double)(tol)==(double)(0) )
2286                {
2287                    tol = eps;
2288                }
2289            }
2290            for(k=2; k<=n; k++)
2291            {
2292                if( iin[k-1]==0 )
2293                {
2294                    y[k] = y[k]-c[k-1]*y[k-1];
2295                }
2296                else
2297                {
2298                    temp = y[k-1];
2299                    y[k-1] = y[k];
2300                    y[k] = temp-c[k-1]*y[k];
2301                }
2302            }
2303            for(k=n; k>=1; k--)
2304            {
2305                if( k<=n-2 )
2306                {
2307                    temp = y[k]-b[k]*y[k+1]-d[k]*y[k+2];
2308                }
2309                else
2310                {
2311                    if( k==n-1 )
2312                    {
2313                        temp = y[k]-b[k]*y[k+1];
2314                    }
2315                    else
2316                    {
2317                        temp = y[k];
2318                    }
2319                }
2320                ak = a[k];
2321                pert = Math.Abs(tol);
2322                if( (double)(ak)<(double)(0) )
2323                {
2324                    pert = -pert;
2325                }
2326                while( true )
2327                {
2328                    absak = Math.Abs(ak);
2329                    if( (double)(absak)<(double)(1) )
2330                    {
2331                        if( (double)(absak)<(double)(sfmin) )
2332                        {
2333                            if( (double)(absak)==(double)(0) | (double)(Math.Abs(temp)*sfmin)>(double)(absak) )
2334                            {
2335                                ak = ak+pert;
2336                                pert = 2*pert;
2337                                continue;
2338                            }
2339                            else
2340                            {
2341                                temp = temp*bignum;
2342                                ak = ak*bignum;
2343                            }
2344                        }
2345                        else
2346                        {
2347                            if( (double)(Math.Abs(temp))>(double)(absak*bignum) )
2348                            {
2349                                ak = ak+pert;
2350                                pert = 2*pert;
2351                                continue;
2352                            }
2353                        }
2354                    }
2355                    break;
2356                }
2357                y[k] = temp/ak;
2358            }
2359        }
2360
2361
2362        private static void internaldlaebz(int ijob,
2363            int nitmax,
2364            int n,
2365            int mmax,
2366            int minp,
2367            double abstol,
2368            double reltol,
2369            double pivmin,
2370            ref double[] d,
2371            ref double[] e,
2372            ref double[] e2,
2373            ref int[] nval,
2374            ref double[,] ab,
2375            ref double[] c,
2376            ref int mout,
2377            ref int[,] nab,
2378            ref double[] work,
2379            ref int[] iwork,
2380            ref int info)
2381        {
2382            int itmp1 = 0;
2383            int itmp2 = 0;
2384            int j = 0;
2385            int ji = 0;
2386            int jit = 0;
2387            int jp = 0;
2388            int kf = 0;
2389            int kfnew = 0;
2390            int kl = 0;
2391            int klnew = 0;
2392            double tmp1 = 0;
2393            double tmp2 = 0;
2394
2395            info = 0;
2396            if( ijob<1 | ijob>3 )
2397            {
2398                info = -1;
2399                return;
2400            }
2401           
2402            //
2403            // Initialize NAB
2404            //
2405            if( ijob==1 )
2406            {
2407               
2408                //
2409                // Compute the number of eigenvalues in the initial intervals.
2410                //
2411                mout = 0;
2412               
2413                //
2414                //DIR$ NOVECTOR
2415                //
2416                for(ji=1; ji<=minp; ji++)
2417                {
2418                    for(jp=1; jp<=2; jp++)
2419                    {
2420                        tmp1 = d[1]-ab[ji,jp];
2421                        if( (double)(Math.Abs(tmp1))<(double)(pivmin) )
2422                        {
2423                            tmp1 = -pivmin;
2424                        }
2425                        nab[ji,jp] = 0;
2426                        if( (double)(tmp1)<=(double)(0) )
2427                        {
2428                            nab[ji,jp] = 1;
2429                        }
2430                        for(j=2; j<=n; j++)
2431                        {
2432                            tmp1 = d[j]-e2[j-1]/tmp1-ab[ji,jp];
2433                            if( (double)(Math.Abs(tmp1))<(double)(pivmin) )
2434                            {
2435                                tmp1 = -pivmin;
2436                            }
2437                            if( (double)(tmp1)<=(double)(0) )
2438                            {
2439                                nab[ji,jp] = nab[ji,jp]+1;
2440                            }
2441                        }
2442                    }
2443                    mout = mout+nab[ji,2]-nab[ji,1];
2444                }
2445                return;
2446            }
2447           
2448            //
2449            // Initialize for loop
2450            //
2451            // KF and KL have the following meaning:
2452            //   Intervals 1,...,KF-1 have converged.
2453            //   Intervals KF,...,KL  still need to be refined.
2454            //
2455            kf = 1;
2456            kl = minp;
2457           
2458            //
2459            // If IJOB=2, initialize C.
2460            // If IJOB=3, use the user-supplied starting point.
2461            //
2462            if( ijob==2 )
2463            {
2464                for(ji=1; ji<=minp; ji++)
2465                {
2466                    c[ji] = 0.5*(ab[ji,1]+ab[ji,2]);
2467                }
2468            }
2469           
2470            //
2471            // Iteration loop
2472            //
2473            for(jit=1; jit<=nitmax; jit++)
2474            {
2475               
2476                //
2477                // Loop over intervals
2478                //
2479                //
2480                // Serial Version of the loop
2481                //
2482                klnew = kl;
2483                for(ji=kf; ji<=kl; ji++)
2484                {
2485                   
2486                    //
2487                    // Compute N(w), the number of eigenvalues less than w
2488                    //
2489                    tmp1 = c[ji];
2490                    tmp2 = d[1]-tmp1;
2491                    itmp1 = 0;
2492                    if( (double)(tmp2)<=(double)(pivmin) )
2493                    {
2494                        itmp1 = 1;
2495                        tmp2 = Math.Min(tmp2, -pivmin);
2496                    }
2497                   
2498                    //
2499                    // A series of compiler directives to defeat vectorization
2500                    // for the next loop
2501                    //
2502                    //*$PL$ CMCHAR=' '
2503                    //CDIR$          NEXTSCALAR
2504                    //C$DIR          SCALAR
2505                    //CDIR$          NEXT SCALAR
2506                    //CVD$L          NOVECTOR
2507                    //CDEC$          NOVECTOR
2508                    //CVD$           NOVECTOR
2509                    //*VDIR          NOVECTOR
2510                    //*VOCL          LOOP,SCALAR
2511                    //CIBM           PREFER SCALAR
2512                    //*$PL$ CMCHAR='*'
2513                    //
2514                    for(j=2; j<=n; j++)
2515                    {
2516                        tmp2 = d[j]-e2[j-1]/tmp2-tmp1;
2517                        if( (double)(tmp2)<=(double)(pivmin) )
2518                        {
2519                            itmp1 = itmp1+1;
2520                            tmp2 = Math.Min(tmp2, -pivmin);
2521                        }
2522                    }
2523                    if( ijob<=2 )
2524                    {
2525                       
2526                        //
2527                        // IJOB=2: Choose all intervals containing eigenvalues.
2528                        //
2529                        // Insure that N(w) is monotone
2530                        //
2531                        itmp1 = Math.Min(nab[ji,2], Math.Max(nab[ji,1], itmp1));
2532                       
2533                        //
2534                        // Update the Queue -- add intervals if both halves
2535                        // contain eigenvalues.
2536                        //
2537                        if( itmp1==nab[ji,2] )
2538                        {
2539                           
2540                            //
2541                            // No eigenvalue in the upper interval:
2542                            // just use the lower interval.
2543                            //
2544                            ab[ji,2] = tmp1;
2545                        }
2546                        else
2547                        {
2548                            if( itmp1==nab[ji,1] )
2549                            {
2550                               
2551                                //
2552                                // No eigenvalue in the lower interval:
2553                                // just use the upper interval.
2554                                //
2555                                ab[ji,1] = tmp1;
2556                            }
2557                            else
2558                            {
2559                                if( klnew<mmax )
2560                                {
2561                                   
2562                                    //
2563                                    // Eigenvalue in both intervals -- add upper to queue.
2564                                    //
2565                                    klnew = klnew+1;
2566                                    ab[klnew,2] = ab[ji,2];
2567                                    nab[klnew,2] = nab[ji,2];
2568                                    ab[klnew,1] = tmp1;
2569                                    nab[klnew,1] = itmp1;
2570                                    ab[ji,2] = tmp1;
2571                                    nab[ji,2] = itmp1;
2572                                }
2573                                else
2574                                {
2575                                    info = mmax+1;
2576                                    return;
2577                                }
2578                            }
2579                        }
2580                    }
2581                    else
2582                    {
2583                       
2584                        //
2585                        // IJOB=3: Binary search.  Keep only the interval
2586                        // containing  w  s.t. N(w) = NVAL
2587                        //
2588                        if( itmp1<=nval[ji] )
2589                        {
2590                            ab[ji,1] = tmp1;
2591                            nab[ji,1] = itmp1;
2592                        }
2593                        if( itmp1>=nval[ji] )
2594                        {
2595                            ab[ji,2] = tmp1;
2596                            nab[ji,2] = itmp1;
2597                        }
2598                    }
2599                }
2600                kl = klnew;
2601               
2602                //
2603                // Check for convergence
2604                //
2605                kfnew = kf;
2606                for(ji=kf; ji<=kl; ji++)
2607                {
2608                    tmp1 = Math.Abs(ab[ji,2]-ab[ji,1]);
2609                    tmp2 = Math.Max(Math.Abs(ab[ji,2]), Math.Abs(ab[ji,1]));
2610                    if( (double)(tmp1)<(double)(Math.Max(abstol, Math.Max(pivmin, reltol*tmp2))) | nab[ji,1]>=nab[ji,2] )
2611                    {
2612                       
2613                        //
2614                        // Converged -- Swap with position KFNEW,
2615                        // then increment KFNEW
2616                        //
2617                        if( ji>kfnew )
2618                        {
2619                            tmp1 = ab[ji,1];
2620                            tmp2 = ab[ji,2];
2621                            itmp1 = nab[ji,1];
2622                            itmp2 = nab[ji,2];
2623                            ab[ji,1] = ab[kfnew,1];
2624                            ab[ji,2] = ab[kfnew,2];
2625                            nab[ji,1] = nab[kfnew,1];
2626                            nab[ji,2] = nab[kfnew,2];
2627                            ab[kfnew,1] = tmp1;
2628                            ab[kfnew,2] = tmp2;
2629                            nab[kfnew,1] = itmp1;
2630                            nab[kfnew,2] = itmp2;
2631                            if( ijob==3 )
2632                            {
2633                                itmp1 = nval[ji];
2634                                nval[ji] = nval[kfnew];
2635                                nval[kfnew] = itmp1;
2636                            }
2637                        }
2638                        kfnew = kfnew+1;
2639                    }
2640                }
2641                kf = kfnew;
2642               
2643                //
2644                // Choose Midpoints
2645                //
2646                for(ji=kf; ji<=kl; ji++)
2647                {
2648                    c[ji] = 0.5*(ab[ji,1]+ab[ji,2]);
2649                }
2650               
2651                //
2652                // If no more intervals to refine, quit.
2653                //
2654                if( kf>kl )
2655                {
2656                    break;
2657                }
2658            }
2659           
2660            //
2661            // Converged
2662            //
2663            info = Math.Max(kl+1-kf, 0);
2664            mout = kl;
2665        }
2666    }
2667}
Note: See TracBrowser for help on using the repository browser.