fd456feeeed5064ae4634671d4b677f74ce6c527
[scilab.git] / scilab / modules / elementary_functions / src / c / expm.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2006 - INRIA - Allan CORNET
4  * Copyright (C) 2012 - Digiteo - Cedric Delamarre
5  *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  *
15  */
16 /*--------------------------------------------------------------------------*/
17
18 #include <stdio.h>
19 #include <string.h>
20 #include "expm.h"
21 #include "basic_functions.h"
22 #include "elementary_functions.h"
23 #include "matrix_left_division.h"
24 #include "matrix_multiplication.h"
25 #include "matrix_right_division.h"
26
27
28 /*
29 purpose
30         compute the exponential of a matrix a by the pade's
31         approximants(subroutine pade).a block diagonalization
32         is performed prior call pade.
33 syntax
34         subroutine dexpm1(ia,n,a,ea,iea,w,iw,ierr)
35
36         integer ia,n,iw,ierr
37         double precision a,ea,w
38         dimension a(ia,n),ea(iea,n),w(*),iw(*)
39
40         ia: the leading dimension of array a.
41         n: the order of the matrices a,ea,x,xi .
42         a: the real double precision array that contains the n*n matrix a
43         ea: the array that contains the n*n exponential of a.
44         iea : the leading dimension of array ea
45         w : work space array of size: n*(2*ia+2*n+5)
46         iw : integer work space array of size 2*n
47         ierr: =0 if the prosessus is normal.
48                 =-1 if n>ia.
49                 =-2 if the block diagonalization is not possible.
50                 =-4 if the subroutine dpade can not computes exp(a)
51 */
52 int dexpms(int _iLeadDim, int _iSize, double *_pdblVal, double *_pdblReturn)
53 {
54     int iErr                    = 0;
55     int iScal                   = 1;
56     int iIndex1                 = 0;
57     int iIndex2                 = 0;
58     int iIndex4                 = 0;
59     int iIndex5                 = 0;
60     int iIndex6                 = 0;
61     int *piWS                   = NULL;
62
63     int iEigenReal              = 0, iEigenImg          = 0;
64     int iBlockStruc             = 0, iRightReduce       = 0;
65     int iInvRightReduce = 0, iScale                     = 0;
66     int iiPadeWS                = 0, idblPadeWS         = 0;
67
68     double *pdblWS              = NULL;
69     double dblANorm             = 0;
70
71     /*global static variable, used to initialization*/
72     sdblExpmN                   = -1;
73
74     for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
75     {
76         double dblAlpha = 0;
77         for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
78         {
79             dblAlpha += dabss(_pdblVal[iIndex2 + _iSize * iIndex1]);
80         }
81         if (dblAlpha > dblANorm)
82         {
83             dblANorm = dblAlpha;
84         }
85     }
86
87     if (dblANorm == 0)
88     {
89         for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
90         {
91             _pdblReturn[iIndex1] = 0;
92             _pdblReturn[iIndex1 + iIndex1 * _iSize] = 1;
93         }
94         return 0;
95     }
96
97     piWS = (int*)malloc(_iSize * _iSize * sizeof(int));
98     pdblWS = (double*)malloc(sizeof(double) * (_iSize * (4 * _iSize + 5)));
99
100     iScale                      = 0;
101     iRightReduce        = _iSize;
102     iInvRightReduce     = iRightReduce + _iLeadDim * _iSize;
103     iEigenReal          = iInvRightReduce + _iLeadDim * _iSize;
104     iEigenImg           = iEigenReal + _iSize;
105     idblPadeWS                  = iEigenImg + _iSize;
106
107     iBlockStruc         = 0;
108     iiPadeWS            = iBlockStruc + _iSize;
109
110     dblANorm = Max(dblANorm, 1);
111     iErr = dbdiaga(_iLeadDim, _iSize, _pdblVal, 0, dblANorm,
112                    &pdblWS[iEigenReal], &pdblWS[iEigenImg], &piWS[iBlockStruc],
113                    &pdblWS[iRightReduce], &pdblWS[iInvRightReduce], &pdblWS[iScale], 1);
114
115     if (iErr)
116     {
117         return -2;
118     }
119
120     for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
121     {
122         vDset(_iSize, 0, &_pdblReturn[iIndex1], _iSize);
123     }
124
125     //compute the pade's approximants of the block.
126     //block origin is shifted before calling pade.
127
128     iIndex1     = 1;
129     iIndex2     = 0;
130
131     //loop on the block.
132     while (iIndex2 <= _iSize)
133     {
134         double dblTemp = 0;
135         double dblAlpha = 0;
136
137         iIndex2 += iIndex1;
138         if (iIndex2 <= _iSize)
139         {
140             iIndex1     = piWS[iBlockStruc + iIndex2 - 1];
141             if (iIndex1 != 1)
142             {
143                 iIndex4 = iIndex2 + iIndex1;
144                 for (iIndex5 = iIndex2 ; iIndex5 < iIndex4 ; iIndex5++)
145                 {
146                     dblTemp += pdblWS[iEigenReal - 1 + iIndex5];
147                 }
148                 dblTemp /= iIndex1;
149
150                 for (iIndex5 = iIndex2 ; iIndex5 < iIndex4 ; iIndex5++)
151                 {
152                     pdblWS[iEigenReal - 1 + iIndex5] -= dblTemp;
153                     _pdblVal[iIndex5 + iIndex5 * _iSize] -= dblTemp;
154                 }
155                 dblAlpha        = 0;
156                 for (iIndex5 = iIndex2 ; iIndex5 < iIndex4 ; iIndex5++)
157                 {
158                     double dblTemp1 = pow(pdblWS[iEigenReal - 1 + iIndex5], 2) + pow(pdblWS[iEigenImg - 1 + iIndex5], 2);
159                     dblTemp1            = dsqrts(dblTemp1);
160                     if (dblTemp1 > dblAlpha)
161                     {
162                         dblAlpha = dblTemp1;
163                     }
164                 }
165
166                 //call pade subroutine.
167                 iErr = dpades(&_pdblVal[iIndex2 + iIndex2 * _iSize], _iSize, iIndex1, &_pdblReturn[iIndex2 + iIndex2 * _iSize],
168                               _iSize, &dblAlpha, &pdblWS[idblPadeWS], &piWS[iiPadeWS]);
169                 /*
170                                                 C2F(pade)(&_pdblVal[iIndex2 + iIndex2 * _iSize], &_iSize, &iIndex1, &_pdblReturn[iIndex2 + iIndex2 * _iSize],
171                                                         &_iSze, &dblAlpha, &pdblWS[idblPadeWS], &piWS[iiPadeWS], &iErr);
172                 */
173                 if (iErr < 0)
174                 {
175                     return 0;
176                 }
177
178                 //remove the effect of origin shift on the block.
179                 dblTemp = dexps(dblTemp);
180                 for (iIndex5 = iIndex2 ; iIndex5 < iIndex4 ; iIndex5++)
181                     for (iIndex6 = iIndex2 ; iIndex6 < iIndex4 ; iIndex6++)
182                     {
183                         _pdblReturn[iIndex5 + iIndex6 * _iSize] *= dblTemp;
184                     }
185             }
186             else
187             {
188                 _pdblReturn[(iIndex2 - 1) + (iIndex2 - 1) * _iSize] =
189                     dexps(_pdblVal[(iIndex2 - 1) + (iIndex2 - 1) * _iSize]);
190             }
191         }
192     }
193
194     //remove the effect of block diagonalization.
195     ddmmuls(&pdblWS[iRightReduce], _iSize, _pdblReturn, _iSize, &pdblWS[idblPadeWS], _iSize, _iSize, _iSize, _iSize);
196     ddmmuls(&pdblWS[idblPadeWS], _iSize, &pdblWS[iInvRightReduce], _iSize, _pdblReturn, _iSize, _iSize, _iSize, _iSize);
197
198     free(piWS);
199     free(pdblWS);
200     return 0;
201 }
202
203 /*
204 purpose
205                 dbdiag reduces a matrix a to block diagonal form by first
206                 reducing it to quasi-triangular form by hqror2 and then by
207                 solving the matrix equation -a11*p+p*a22=a12 to introduce zeros
208                 above the diagonal.
209                 right transformation is factored : p*d*u*y ;where:
210                         p is a permutation matrix and d positive diagonal matrix
211                         u is orthogonal and y block upper triangular with identity
212                                 blocks on the diagonal
213
214         starred parameters are altered by the subroutine
215
216
217                 *a              an array that initially contains the m x n matrix
218                                 to be reduced. on return,  see job
219
220                 lda             the leading dimension of array a. and array x,xi.
221                 n               the order of the matrices a,x,xi
222
223                 epsshr  the minimal conditionnement allowed for linear sytems
224
225                 rmax    the maximum size allowed for any element of the
226                                 transformations.
227
228                 *er             a singly subscripted real array containing the real
229                                 parts of the eigenvalues.
230
231                 *ei             a singly subscripted real array containing the imaginary
232                                 parts of the eigenvalues.
233
234                 *bs             a singly subscripted integer array that contains block
235                                 structure information.  if there is a block of order
236                                 k starting at a(l,l) in the output matrix a, then
237                                 bs(l) contains the positive integer k, bs(l+1) contains
238                                 -(k-1), bs(la+2) = -(k-2), ..., bs(l+k-1) = -1.
239                                 thus a positive integer in the l-th entry of bs
240                                 indicates a new block of order bs(l) starting at a(l,l).
241
242                 *x              contains,  either right reducing transformation u*y,
243                                 either orthogonal tranformation u (see job)
244
245                 *xi             xi contains the inverse reducing matrix transformation
246                                 or y matrix (see job)
247
248                 *scale  contains the scale factor and definitions of p size(n)
249
250                 job             integer parametre specifying outputed transformations
251                                 job=0 : a contains block diagonal form
252                                                 x right transformation
253                                                 xi dummy variable
254                                 job=1 : like job=0 and xi contain x**-1
255                                 job=2   a contains block diagonal form
256                                                 x contains u and xi contain y
257                                 job=3 : a contains:
258                                                         -block diagonal form in the diagonal blocks
259                                                         -a factorisation of y in the upper triangular
260                                                 x contains u
261                                                 xi dummy
262                                 *fail   a logical variable which is false on normal return and
263                                                 true if there is any error in bdiag.
264 */
265 int dbdiaga(int _iLeadDim, int _iSize, double *_pdblVal, double _dblEps,
266             double _dblMax, double *_pdblEigenReal, double *_pdblEigenImg,
267             int *_piBlockStruc, double *_pdblRightReduce, double *_pdblInvRightReduce,
268             double *_pdblScale, int _iMode)
269 {
270     int iErr                    = 0;
271     int iIndex1                 = 0;
272     int iIndex2                 = 0;
273     int iIndex3                 = 0;
274     int iLow                    = 0;
275     int iHigh                   = 0;
276     int iLoop11                 = 0;
277     int iLoop22                 = 0;
278     int iSize11                 = 0;
279     int iSize22                 = 0;
280     int iLoop22m1               = 0;
281     int iK                              = 0;
282     int iL                              = 0;
283     int inK                             = 0;
284     int iKm1                    = 0;
285     int iKm2                    = 0;
286
287     int iMinusOne               = -1;
288     int iZero                   = 0;
289     int iOne                    = 1;
290     int iTwo                    = 2;
291
292     double dblMinusOne  = -1;
293     double dblZero              = 0;
294     double dblOne               = 1;
295     double dblTwo               = 2;
296
297     double dblEps               = 0;
298     double dblAvgReal   = 0;
299     double dblAvgImg    = 0;
300     double dblD                 = 0;
301
302     for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
303     {
304         double dblTemp  = 0;
305         for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
306         {
307             dblTemp += dabss(_pdblVal[iIndex2 + _iSize * iIndex1]);
308         }
309
310         dblEps = Max(dblEps, dblTemp);
311     }
312     if (dblEps == 0)
313     {
314         dblEps = 1;
315     }
316
317     dblEps = nc_eps_machine() * dblEps;
318
319     //convert a to upper hessenberg form.
320     dbalancs(_iLeadDim, _iSize, _pdblVal, &iLow, &iHigh, _pdblScale);
321     dorthess(_iLeadDim, _iSize, iLow, iHigh, _pdblVal, _pdblEigenReal);
322     dortrans(_iLeadDim, _iSize, iLow, iHigh, _pdblVal, _pdblEigenReal, _pdblRightReduce);
323
324     //convert a to quasi-upper triangular form by qr method.
325     iErr = dhqror2s(_iLeadDim, _iSize, 1, _iSize, _pdblVal, _pdblEigenReal, _pdblEigenImg, _pdblRightReduce, 21);
326
327     //check to see if hqror2 failed in computing any eigenvalue
328     if (iErr != 0)
329     {
330         return iErr;
331     }
332
333     /*
334         reduce a to block diagonal form
335
336
337         segment a into 4 matrices: a11, a da11 x da11 block
338         whose (1,1)-element is at a(l11,l11))  a22, a da22 x da22
339         block whose (1,1)-element is at a(l22,l22)) a12,
340         a da11 x da22 block whose (1,1)-element is at a(l11,l22))
341         and a21, a da22 x da11 block = 0 whose (1,1)-
342         element is at a(l22,l11).
343
344
345
346         this loop uses l11 as loop index and splits off a block
347         starting at a(l11,l11).
348     */
349     iLoop11     = 1;
350
351     //L40
352     while (iLoop11 <= _iSize)
353     {
354         iLoop22 = iLoop11;
355
356         //this loop uses da11 as loop variable and attempts to split
357         //off a block of size da11 starting at a(l11,l11)
358
359         //sorry for "while(1)" but the original source code in fortran
360         //do loop until goto outside the loop
361         //L50
362         while (1)
363         {
364             if (iLoop22 == iLoop11)
365             {
366                 iSize11         = 1;
367                 if (iLoop11 != (_iSize - 1))
368                     if (dabss(_pdblVal[iLoop11 + (iLoop11 - 1) * _iSize]) > dblEps)
369                     {
370                         iSize11 = 2;
371                     }
372
373                 iLoop22         = iLoop11 + iSize11;
374                 iLoop22m1       = iLoop22 - 1;
375
376             }
377             else
378             {
379                 //compute the average of the eigenvalues in a11
380                 dblAvgReal      = 0;
381                 dblAvgImg       = 0;
382                 for (iIndex1 = iLoop11 - 1 ; iIndex1 < iLoop22m1 - 1 ; iIndex1++)
383                 {
384                     dblAvgReal  += _pdblEigenReal[iIndex1];
385                     dblAvgImg   += dabss(_pdblEigenImg[iIndex1]);
386                 }
387                 dblAvgReal      /= iSize11;
388                 dblAvgImg       /= iSize11;
389
390                 //loop on eigenvalues of a22 to find the one closest to the av
391                 dblD            = pow(dblAvgReal - _pdblEigenReal[iLoop22], 2) + pow(dblAvgImg - _pdblEigenImg[iLoop22], 2);
392                 iK                      = iLoop22;
393                 iL                      = iK + 1;
394
395                 if (iLoop22 != (_iSize - 1))
396                     if (dabss(_pdblVal[(iLoop22 + 1) + iLoop22 * _iSize]) > dblEps)
397                     {
398                         iL      = iLoop22 + 2;
399                     }
400
401                 //L71
402
403                 //L80
404                 while (iL < _iSize)
405                 {
406                     double dblTemp = pow(dblAvgReal - _pdblEigenReal[iL], 2) + pow(dblAvgImg - _pdblEigenImg[iL], 2);
407                     if (dblTemp < dblD)
408                     {
409                         iK              = iL;
410                         dblD    = dblTemp;
411                     }
412                     iL++;
413                     if (iL <= _iSize)
414                         if (dabss(_pdblVal[iL + (iL * - 1) * _iSize]) > dblEps)
415                         {
416                             iL++;
417                         }
418                 }
419                 //loop to move the eigenvalue just located
420                 //into first position of block a22.
421                 if (iK != (_iSize - 1) && dabss(_pdblVal[(iK + 1) + iK * _iSize] > dblEps))
422                 {
423                     inK = 1;
424                     while (iK > iLoop22)
425                     {
426                         iKm1 = iK - 1;
427                         if (dabss(_pdblVal[iKm1 + (iK - 2) * _iSize]) >= dblEps)
428                         {
429                             double dblE1 = 0, dblE2 = 0;
430                             //we're swapping the closest block with a 2 x 2
431                             iKm2 = iK - 2;
432                             //C2F(exch)(&_iLeadDim, &_iSize, _pdblVal, _pdblRightReduce, &iKm2, &iTwo, &iOne);
433                             dexchs(_iLeadDim, _iSize, _pdblVal, _pdblRightReduce, iKm2, 2, 1);
434                             //try to split this block into 2 real eigenvalues
435                             //C2F(split)(_pdblVal, _pdblRightReduce, &_iSize, &iKm1, &dblE1, &dblE2, &_iLeadDim, &_iLeadDim);
436                             dsplits(_pdblVal, _pdblRightReduce, _iSize, iKm1, &dblE1, &dblE2, _iLeadDim, _iLeadDim);
437                         }
438                     }
439                 }
440
441                 inK     = 2;
442             }
443
444             //L240
445             if (iLoop22 <= _iSize)
446             {
447                 //attempt to split off a block of size da11.
448                 iSize22 = _iSize - iLoop22 + 1;
449                 //save a12 in its transpose form in block a21.
450
451                 for (iIndex1 = iLoop11 - 1 ; iIndex1 < iLoop22m1 ; iIndex1++)
452                 {
453                     for (iIndex2 = iLoop22 - 1 ; iIndex2 < _iSize ; iIndex2++)
454                     {
455                         _pdblVal[iIndex2 + iIndex1 * _iSize] = _pdblVal[iIndex1 + iIndex2 * _iSize];
456                     }
457                 }
458
459                 //convert a11 to lower quasi-triangular and multiply it by -1 and
460                 //a12 appropriately (for solving -a11*p+p*a22=a12).
461                 C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &_iSize, &dblOne, &iZero);
462                 C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &iLoop22m1, &dblMinusOne, &iOne);
463
464                 //solve -a11*p + p*a22 = a12.
465
466                 C2F(shrslv)(&_pdblVal[(iLoop11 - 1) + (iLoop11 - 1) * _iSize], &_pdblVal[(iLoop22 - 1) + (iLoop22 - 1) * _iSize],
467                             &_pdblVal[(iLoop11 - 1) + (iLoop22 - 1) * _iSize], &iSize11, &iSize22,
468                             &_iLeadDim, &_iLeadDim, &_iLeadDim, &dblEps, &_dblEps, &_dblMax, &iErr);
469
470                 if (iErr)
471                 {
472                     //change a11 back to upper quasi-triangular.
473                     C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &iLoop22m1, &dblOne, &iOne);
474                     C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &iLoop22m1, &dblMinusOne, &iOne);
475
476                     //was unable to solve for p - try again
477                     //move saved a12 back into its correct position.
478                     for (iIndex1 = iLoop11 ; iIndex1 < iLoop22m1 ; iIndex1++)
479                     {
480                         for (iIndex2 = iLoop22 ; iIndex2 < _iSize ; iIndex2++)
481                         {
482                             _pdblVal[iIndex1 + iIndex2 * _iSize]        = _pdblVal[iIndex2 + iIndex1 * _iSize];
483                             _pdblVal[iIndex2 + iIndex1 * _iSize]        = 0;
484                         }
485                     }
486                 }//if iErr
487                 else
488                 {
489                     break;
490                 }
491             }//if iLoop22 <= _iSize
492             else
493             {
494                 break;
495             }
496         }// while 1
497         //L290
498         //change solution to p to proper form.
499         if (iLoop22 <= _iSize)
500         {
501             C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &_iSize, &dblOne, &iZero);
502             C2F(dad)(_pdblVal, &_iLeadDim, &iLoop11, &iLoop22m1, &iLoop11, &iLoop22m1, &dblMinusOne, &iOne);
503         }
504
505         //store block size in array bs.
506         _piBlockStruc[iLoop11 - 1]      = iSize11;
507         iIndex1 = iSize11 - 1;
508         if (iIndex1 != 0)
509         {
510             for (iIndex2 = 0 ; iIndex2 < iIndex1; iIndex2++)
511             {
512                 _piBlockStruc[(iLoop11 - 1) + iIndex2] = -(iSize11 - iIndex2);
513             }
514         }
515         iLoop11 = iLoop22;
516     }//while iLoop11 < _iSize (L350)
517
518     //set transformations matrices as required
519     if (_iMode == 3)
520     {
521         return 0;
522     }
523
524     switch (_iMode)
525     {
526         case 2 :        //extract non orthogonal tranformation from a
527             for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
528             {
529                 vDset(_iSize, 0, &_pdblInvRightReduce[iIndex1 * _iSize], 1);
530             }
531             vDset(_iSize, 1, _pdblInvRightReduce, _iLeadDim + 1);
532
533             iLoop22 = 1;
534             while (iLoop11 <= _iSize)
535             {
536                 iLoop11 = iLoop22;
537                 iLoop22 = iLoop11 + _piBlockStruc[iLoop11];
538                 for (iIndex1 = (iLoop22 - 1) ; iIndex1 < _iSize ; iIndex1++)
539                 {
540                     for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
541                     {
542                         int iTemp1 = iLoop22 - iLoop11;
543                         _pdblInvRightReduce[iIndex2 + iIndex1 * _iSize] +=
544                             C2F(ddot)(&iTemp1, &_pdblInvRightReduce[iIndex2 + (iLoop11 - 1) * _iSize], &_iLeadDim, &_pdblVal[(iLoop11 - 1) + iIndex1 * _iSize], &iOne);
545                     }
546                 }
547             }
548             break;
549         case 1 :        //compute inverse transformation
550             for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
551             {
552                 for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
553                 {
554                     _pdblInvRightReduce[iIndex1 + iIndex2 * _iSize] = _pdblRightReduce[iIndex2 + iIndex1 * _iSize];
555                 }
556             }
557             iLoop22 = 1;
558
559             while (iLoop22 <= _iSize)
560             {
561                 iLoop11 = iLoop22;
562                 iLoop22 = iLoop11 + _piBlockStruc[iLoop11 - 1];
563                 if (iLoop22 <= _iSize)
564                 {
565                     iLoop22m1   = iLoop22 - 1;
566
567                     for (iIndex1 = (iLoop11 - 1) ; iIndex1 < iLoop22m1 ; iIndex1++)
568                     {
569                         for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
570                         {
571                             int iTemp1 = _iSize - iLoop22m1;
572                             _pdblInvRightReduce[iIndex1 + iIndex2 * _iSize] -=
573                                 C2F(ddot)(&iTemp1, &_pdblVal[iIndex1 + (iLoop22 - 1) * _iSize], &_iLeadDim, &_pdblInvRightReduce[(iLoop22 - 1) + iIndex2 * _iSize], &iOne);
574                         }
575                     }
576                 }
577             }//while(iLoop22 <= _iSize)
578
579             if (iHigh != iLow)
580             {
581                 for (iIndex1 = iLow ; iIndex1 < iHigh ; iIndex1++)
582                 {
583                     double dblTemp = 1.0 / _pdblScale[iIndex1];
584                     for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
585                     {
586                         _pdblInvRightReduce[iIndex2 + iIndex1 * _iSize] *= dblTemp;
587                     }
588                 }
589             }
590
591             for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
592             {
593                 iIndex2 = iIndex1;
594                 if (iIndex2 < iLow || iIndex2 > iHigh)
595                 {
596                     if (iIndex2 < iLow)
597                     {
598                         iIndex2 = iLow - iIndex1;
599                     }
600
601                     iK  = (int)_pdblScale[iIndex2];
602                     if (iK != iIndex2)
603                     {
604                         for (iIndex3 = 0 ; iIndex3 < _iSize ; iIndex3++)
605                         {
606                             vSwitchVal(_pdblInvRightReduce, iIndex3 + iIndex2 * _iSize, iIndex3 + iK * _iSize);
607                         }
608                     }
609                 }
610             }
611         // Warning, case 1 continue to default process
612         //break;
613         default :       //compute right transformation
614             iLoop22     = 1;
615             while (iLoop22 <= _iSize)
616             {
617                 iLoop11 = iLoop22;
618                 iLoop22 = iLoop11 + _piBlockStruc[iLoop11 - 1];
619                 if (iLoop11 <= _iSize)
620                 {
621                     for (iIndex1 = (iLoop22 - 1) ; iIndex1 < _iSize ; iIndex1++)
622                     {
623                         for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
624                         {
625                             int iTemp1 = iLoop22 - iLoop11;
626                             _pdblRightReduce[iIndex2 + iIndex1 * _iSize] +=
627                                 C2F(ddot)(&iTemp1, &_pdblRightReduce[iIndex2 + (iLoop11 - 1) * _iSize], &_iLeadDim, &_pdblVal[(iLoop11 - 1) + iIndex1 * _iSize], &iOne);
628                         }
629                     }
630                 }
631             }
632
633             C2F(balbak)(&_iLeadDim, &_iSize, &iLow, &iHigh, _pdblScale, &_iSize, _pdblRightReduce);
634             break;
635     }
636
637     //set zeros in the matrix a
638     iLoop11     = 1;
639     iLoop22 = iLoop11 + _piBlockStruc[iLoop11 - 1];
640     while (iLoop22 < _iSize)
641     {
642         iLoop22 = iLoop11 + _piBlockStruc[iLoop11 - 1];
643         if (iLoop22 <= _iSize)
644         {
645             iLoop22m1   = iLoop22 - 1;
646             for (iIndex1 = (iLoop11 - 1) ; iIndex1 < iLoop22m1 ; iIndex1++)
647             {
648                 vDset(_iSize - iLoop22m1, 0, &_pdblVal[iIndex1 + (iLoop22 - 1) * _iSize], _iLeadDim);
649                 vDset(_iSize - iLoop22m1, 0, &_pdblVal[(iLoop22 - 1) + iIndex1 * _iSize], 1);
650             }
651             iLoop11 = iLoop22;
652         }
653     }
654
655
656     return iErr;
657 }
658
659 /*
660 purpose
661
662      this subroutine balances a real matrix and isolates
663      eigenvalues whenever possible.
664 syntax
665
666      on input:
667
668         nm must be set to the row dimension of two-dimensional
669           array parameters as declared in the calling program
670           dimension statement;
671
672         n is the order of the matrix;
673
674         a contains the input matrix to be balanced.
675
676      on output:
677
678         a contains the balanced matrix;
679
680         low and igh are two integers such that a(i,j)
681           is equal to zero if
682            (1) i is greater than j and
683            (2) j=1,...,low-1 or i=igh+1,...,n;
684
685         scale contains information determining the
686            permutations and scaling factors used.
687
688      suppose that the principal submatrix in rows low through igh
689      has been balanced, that p(j) denotes the index interchanged
690      with j during the permutation step, and that the elements
691      of the diagonal matrix used are denoted by d(i,j).  then
692         scale(j) = p(j),    for j = 1,...,low-1
693                  = d(j,j),      j = low,...,igh
694                  = p(j)         j = igh+1,...,n.
695      the order in which the interchanges are made is n to igh+1,
696      then 1 to low-1.
697
698      note that 1 is returned for igh if igh is zero formally.
699
700      the algol procedure exc contained in balance appears in
701      balanc  in line.  (note that the algol roles of identifiers
702      k,l have been reversed.)
703
704 originator
705
706      this subroutine is a translation of the algol procedure balance,
707      num. math. 13, 293-304(1969) by parlett and reinsch.
708      handbook for auto. comp., vol.ii-linear algebra, 315-326(1971).
709
710      questions and comments should be directed to b. s. garbow,
711      applied mathematics division, argonne national laboratory
712
713      ------------------------------------------------------------------
714
715      :::::::::: radix is a machine dependent parameter specifying
716                 the base of the machine floating point representation.
717 */
718 int dbalancs(int _iRows, int _iSize, double *_pdblVal, int *_piLow, int *_piHigh, double *_pdblScale)
719 {
720     double dblRadix             = 2;
721     double dblRadix2    = dblRadix * dblRadix;
722     int iIndex1                 = 0;
723     int iIndex2                 = 0;
724     int iLoop1                  = _iSize;
725     int iLoop2                  = 1;
726
727     //flag
728     int bLoop1                  = 1;
729     int bLoop2                  = 1;
730     int iNoConv                 = 0;
731
732     double dblRowSum    = 0;
733     double dblColSum    = 0;
734     double dblG                 = 0;
735     double dblF                 = 0;
736     double dblS                 = 0;
737
738     //go to 100
739     while (bLoop1)
740     {
741         int iTemp = 0;
742         bLoop1 = 0;
743         for (iIndex1 = 0 ; iIndex1 < iLoop1 ; iIndex1++)
744         {
745             iTemp = iLoop1 - iIndex1 - 1;
746             for (iIndex2 = iLoop2 - 1 ; iIndex2 < iLoop1 ; iIndex2++)
747             {
748                 if (iIndex2 == iTemp)
749                 {
750                     continue;
751                 }
752                 if (_pdblVal[iTemp + iIndex2 * _iSize] != 0)
753                 {
754                     break;
755                 }
756             }
757             if (iIndex2 < iLoop1 && iIndex1 < iLoop1 && _pdblVal[iTemp + iIndex2 * _iSize] != 0)
758             {
759                 continue;
760             }
761             bLoop1 = 1;
762             break;
763         }
764         if (bLoop1)
765         {
766             vExchangeVal(       _pdblScale, _pdblVal,
767                             0, iLoop1 - 1,
768                             iLoop2 - 1, _iSize,
769                             _iSize, iTemp, iLoop1 - 1);
770             if (iLoop1 == 1)
771             {
772                 *_piLow         = iLoop2;
773                 *_piHigh        = iLoop1;
774                 return 0;
775             }
776             else
777             {
778                 iLoop1--;
779             }
780         }
781     }
782
783     while (bLoop2)
784     {
785         int iTemp = 0;
786         bLoop2 = 0;
787         for (iIndex1 = iLoop2 - 1 ; iIndex1 < iLoop1 ; iIndex1++)
788         {
789             for (iIndex2 = iLoop2 - 1 ; iIndex2 < iLoop1 ; iIndex2++)
790             {
791                 if (iIndex2 == iIndex1)
792                 {
793                     continue;
794                 }
795                 if (_pdblVal[iIndex2 + iIndex1 * _iSize] != 0)
796                 {
797                     break;
798                 }
799             }
800             if (iIndex2 < iLoop1 && iIndex1 < iLoop1 && _pdblVal[iIndex2 + iIndex1 * _iSize] != 0)
801             {
802                 continue;
803             }
804             bLoop2 = 1;
805             break;
806         }
807         if (bLoop2)
808         {
809             vExchangeVal(       _pdblScale, _pdblVal,
810                             0, iLoop1 - 1,
811                             iLoop2 - 1, _iSize,
812                             _iSize, iIndex1, iLoop2 - 1);
813
814             iLoop2++;
815         }
816     }
817
818     for (iIndex1 = iLoop2 - 1 ; iIndex1 < iLoop1 ; iIndex1++)
819     {
820         _pdblScale[iIndex1] = 1;
821     }
822
823     do
824     {
825         iNoConv = 0;
826         for (iIndex1 = iLoop2 - 1 ; iIndex1 < iLoop1 ; iIndex1++)
827         {
828             dblRowSum = 0;
829             dblColSum = 0;
830
831             for (iIndex2 = iLoop2 - 1 ; iIndex2 < iLoop1 ; iIndex2++)
832             {
833                 dblColSum += dabss(_pdblVal[iIndex2 + iIndex1 * _iSize]);
834                 dblRowSum += dabss(_pdblVal[iIndex1 + iIndex2 * _iSize]);
835             }
836
837             //:::::::::: guard against zero c or r due to underflow ::::::::::
838             if (dblRowSum == 0 || dblColSum == 0)
839             {
840                 continue;
841             }
842
843             dblG        = dblRowSum / dblRadix;
844             dblF        = 1;
845             dblS        = dblColSum + dblRowSum;
846
847             while (dblColSum < dblG)
848             {
849                 dblF            *= dblRadix;
850                 dblColSum       *= dblRadix2;
851             }
852
853             dblG        = dblRowSum * dblRadix;
854
855             while (dblColSum >= dblG)
856             {
857                 dblF            /= dblRadix;
858                 dblColSum       /= dblRadix2;
859             }
860
861             //:::::::::: now balance ::::::::::
862             if ( ((dblColSum + dblRowSum) / dblF) >= (0.95 * dblS))
863             {
864                 continue;
865             }
866
867             dblG                                = 1 / dblF;
868             _pdblScale[iIndex1] *= dblF;
869             iNoConv                             = 1;
870
871             for (iIndex2 = iLoop2 - 1 ; iIndex2 < _iSize ; iIndex2++)
872             {
873                 _pdblVal[iIndex1 + iIndex2 * _iSize]    *= dblG;
874             }
875
876             for (iIndex2 = 0 ; iIndex2 < iLoop1 ; iIndex2++)
877             {
878                 _pdblVal[iIndex2 + iIndex1 * _iSize]    *= dblF;
879             }
880         }
881     }
882     while (iNoConv == 1);
883
884     *_piLow             = iLoop2; //Check if this value is use like array index or not ( 0 indexed or 1 indexed )
885     *_piHigh    = iLoop1; //Check if this value is use like array index or not ( 0 indexed or 1 indexed )
886
887     //:::::::::: last card of balanc ::::::::::
888     return 0;
889 }
890
891 void vExchangeVal(double *_pdblScale, double *_pdblVal,
892                   int _iStart1, int _iEnd1,
893                   int _iStart2, int _iEnd2,
894                   int _iSize, int _iCoord1, int _iCoord2)
895 {
896     int iIndex                  = 0;
897
898     _pdblScale[_iEnd1]  = _iCoord1;
899     if (_iCoord1 == _iEnd1)
900     {
901         return;
902     }
903     for (iIndex = _iStart1 ; iIndex <= _iEnd1 ; iIndex++)
904         vSwitchVal(             _pdblVal,
905                         iIndex + _iCoord1 * _iSize,
906                         iIndex + _iCoord2 * _iSize);
907
908     for (iIndex = _iStart2 ; iIndex < _iEnd2 ; iIndex++)
909         vSwitchVal(             _pdblVal,
910                         _iCoord1 + iIndex * _iSize,
911                         _iCoord2 + iIndex * _iSize);
912 }
913
914 void vSwitchVal(double *_pdblVal, int _iPos1, int _iPos2)
915 {
916     double dblTemp              = 0;
917
918     dblTemp                             = _pdblVal[_iPos1];
919     _pdblVal[_iPos1]    = _pdblVal[_iPos2];
920     _pdblVal[_iPos2]    = dblTemp;
921 }
922
923 /*
924  purpose
925
926     given a real general matrix, this subroutine
927     reduces a submatrix situated in rows and columns
928     low through igh to upper hessenberg form by
929     orthogonal similarity transformations.
930
931  syntax
932
933      subroutine orthes(nm,n,low,igh,a,ort)
934
935     on input:
936
937        nm must be set to the row dimension of two-dimensional
938          array parameters as declared in the calling program
939          dimension statement;
940
941        n is the order of the matrix;
942
943        low and igh are integers determined by the balancing
944          subroutine  balanc.  if  balanc  has not been used,
945          set low=1, igh=n;
946
947        a contains the input matrix.
948
949     on output:
950
951        a contains the hessenberg matrix.  information about
952          the orthogonal transformations used in the reduction
953          is stored in the remaining triangle under the
954          hessenberg matrix;
955
956        ort contains further information about the transformations.
957          only elements low through igh are used.
958
959 originator
960
961     this subroutine is a translation of the algol procedure orthes,
962     num. math. 12, 349-368(1968) by martin and wilkinson.
963     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
964     questions and comments should be directed to b. s. garbow,
965     applied mathematics division, argonne national laboratory
966 */
967 int dorthess(int _iLead, int _iSize, int _iLow, int _iHigh, double *_pdblVal, double *_pdblOrt)
968 {
969     int iIndex1         = 0;
970     int iIndex2         = 0;
971     int iIndex3         = 0;
972
973     int iNewHigh        = _iHigh - 1;
974     int iNewLow         = _iLow + 1;
975
976     if (iNewHigh < iNewLow)
977     {
978         return 0;
979     }
980
981     for (iIndex1 = iNewLow - 1 ; iIndex1 < iNewHigh ; iIndex1++)
982     {
983         double dblH                     = 0;
984         double dblG                     = 0;
985         double dblScale         = 0;
986         int iOffMax                     = 0;
987
988         _pdblOrt[iIndex1]       = 0;
989
990         //:::::::::: scale column (algol tol then not needed) ::::::::::
991         for (iIndex2 = iIndex1 ; iIndex2 < _iHigh; iIndex2++)
992         {
993             dblScale += dabss(_pdblVal[iIndex2 + (iIndex1 - 1) * _iSize]);
994         }
995
996         if (dblScale == 0)
997         {
998             continue;
999         }
1000
1001         iOffMax = iIndex1 + _iHigh;
1002
1003         //:::::::::: for i=igh step -1 until m do -- ::::::::::
1004         for (iIndex2 = iIndex1 ; iIndex2 < _iHigh; iIndex2++)
1005         {
1006             iIndex3                             = iOffMax - (iIndex2 + 1);
1007             _pdblOrt[iIndex3]   = _pdblVal[iIndex3 + (iIndex1 - 1) * _iSize] / dblScale;
1008             dblH                                += _pdblOrt[iIndex3] * _pdblOrt[iIndex3];
1009         }
1010
1011         dblG                            = -dsigns(dsqrts(dblH), _pdblOrt[iIndex1]);
1012         dblH                            -= _pdblOrt[iIndex1] * dblG;
1013         _pdblOrt[iIndex1]       -= dblG;
1014
1015         //:::::::::: form (i-(u*ut)/h) * a ::::::::::
1016         for (iIndex2 = iIndex1 ; iIndex2 < _iSize; iIndex2++)
1017         {
1018             double dblF = 0;
1019
1020             //:::::::::: for i=igh step -1 until m do -- ::::::::::
1021             for (iIndex3 = iIndex1; iIndex3 < _iHigh; iIndex3++)
1022             {
1023                 int iTemp       = iOffMax - (iIndex3 + 1);
1024                 dblF            += _pdblOrt[iTemp] * _pdblVal[iTemp + iIndex2 * _iSize];
1025             }
1026
1027             dblF /= dblH;
1028
1029             for (iIndex3 = iIndex1; iIndex3 < _iHigh; iIndex3++)
1030             {
1031                 _pdblVal[iIndex3 + iIndex2 * _iSize] -= dblF * _pdblOrt[iIndex3];
1032             }
1033         }
1034
1035         //:::::::::: form (i-(u*ut)/h)*a*(i-(u*ut)/h) ::::::::::
1036         for (iIndex2 = 0 ; iIndex2 < _iHigh; iIndex2++)
1037         {
1038             double dblF = 0;
1039
1040             //:::::::::: for j=igh step -1 until m do -- ::::::::::
1041             for (iIndex3 = iIndex1; iIndex3 < _iHigh; iIndex3++)
1042             {
1043                 int iTemp       = iOffMax - (iIndex3 + 1);
1044                 dblF            += _pdblOrt[iTemp] * _pdblVal[iIndex2 + iTemp * _iSize];
1045             }
1046
1047             dblF /= dblH;
1048
1049             for (iIndex3 = iIndex1; iIndex3 < _iHigh; iIndex3++)
1050             {
1051                 _pdblVal[iIndex2 + iIndex3 * _iSize] -= dblF * _pdblOrt[iIndex3];
1052             }
1053         }
1054
1055         _pdblOrt[iIndex1]                                                       *= dblScale;
1056         _pdblVal[iIndex1 + (iIndex1 - 1) * _iSize]      = dblScale * dblG;
1057     }
1058     //:::::::::: last card of orthes ::::::::::
1059     return 0;
1060 }
1061
1062 /*
1063 purpose
1064
1065     this subroutine accumulates the orthogonal similarity
1066     transformations used in the reduction of a real general
1067     matrix to upper hessenberg form by  orthes.
1068
1069 syntax
1070
1071      subroutine ortran(nm,n,low,igh,a,ort,z)
1072
1073     on input:
1074
1075        nm must be set to the row dimension of two-dimensional
1076          array parameters as declared in the calling program
1077          dimension statement;
1078
1079        n is the order of the matrix;
1080
1081        low and igh are integers determined by the balancing
1082          subroutine  balanc.  if  balanc  has not been used,
1083          set low=1, igh=n;
1084
1085        a contains information about the orthogonal trans-
1086          formations used in the reduction by  orthes
1087          in its strict lower triangle;
1088
1089        ort contains further information about the trans-
1090          formations used in the reduction by  orthes.
1091          only elements low through igh are used.
1092
1093     on output:
1094
1095        z contains the transformation matrix produced in the
1096          reduction by  orthes;
1097
1098        ort has been altered.
1099
1100 originator
1101
1102     this subroutine is a translation of the algol procedure ortrans,
1103     num. math. 16, 181-204(1970) by peters and wilkinson.
1104     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
1105
1106     questions and comments should be directed to b. s. garbow,
1107     applied mathematics division, argonne national laboratory
1108 */
1109 int dortrans(int _iLead, int _iSize, int _iLow, int _iHigh, double *_pdblVal, double *_pdblOrt, double *_pdblTrans)
1110 {
1111     int iIndex1 = 0;
1112     int iIndex2 = 0;
1113     int iIndex3 = 0;
1114     int iLoop = 0;
1115
1116     for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
1117     {
1118         for (iIndex2 = 0 ; iIndex2 < _iSize ; iIndex2++)
1119         {
1120             _pdblTrans[iIndex1 + iIndex2 * _iSize] = 0;
1121         }
1122
1123         _pdblTrans[iIndex1 + iIndex1 * _iSize] = 1;
1124     }
1125
1126     iLoop = _iHigh - _iLow - 1;
1127     if (iLoop < 0)
1128     {
1129         return 0;
1130     }
1131
1132     //:::::::::: for mp=igh-1 step -1 until low+1 do -- ::::::::::
1133     for (iIndex1 = 0 ; iIndex1 < iLoop ; iIndex1++)
1134     {
1135         int iTemp       = _iHigh - (iIndex1 + 1) - 1;
1136
1137         if (_pdblVal[iTemp + (iTemp - 1) * _iSize] == 0)
1138         {
1139             continue;
1140         }
1141
1142         for (iIndex2 = iTemp + 1 ;  iIndex2 < _iHigh ; iIndex2++)
1143         {
1144             _pdblOrt[iIndex2] = _pdblVal[iIndex2 + (iTemp - 1) * _iSize];
1145         }
1146
1147         for (iIndex2 = iTemp ;  iIndex2 < _iHigh ; iIndex2++)
1148         {
1149             double dblG = 0;
1150
1151             for (iIndex3 = iTemp ;  iIndex3 < _iHigh ; iIndex3++)
1152             {
1153                 dblG    += _pdblOrt[iIndex3] * _pdblTrans[iIndex3 + iIndex2 * _iSize];
1154             }
1155
1156             //::::::::::divisor below is negative of h formed in orthes.
1157             //::::::::::double division avoids possible underflow
1158             dblG        = (dblG / _pdblOrt[iTemp] ) / _pdblVal[iTemp + (iTemp - 1) * _iSize];
1159
1160             for (iIndex3 = iTemp ;  iIndex3 < _iHigh ; iIndex3++)
1161             {
1162                 _pdblTrans[iIndex3 + iIndex2 * _iSize]  += dblG * _pdblOrt[iIndex3];
1163             }
1164
1165         }
1166     }
1167
1168     return 0;
1169 }
1170
1171 /*
1172
1173         this subroutine is a translation of the algol procedure hqr2,
1174         num. math. 16, 181-204(1970) by peters and wilkinson.
1175         handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
1176
1177 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1178
1179     MODIFICATIONS WRT EISPACK VERSION
1180     ---------------------------------
1181       1. 1x1 and 2x2 diagonal blocks are clearly isolated by
1182          forcing subdiagonal entries to zero
1183       2. Merging of hqr/hqr2 driven by a job parameter
1184
1185 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1186
1187     This subroutine finds the eigenvalues of a real upper
1188     hessenberg matrix by the qr method. In addition, the
1189     orthogonal transformation leading to the Schur form is
1190     accumulated
1191
1192     on input
1193
1194        nm must be set to the row dimension of two-dimensional
1195          array parameters as declared in the calling program
1196          dimension statement.
1197
1198        n is the order of the matrix.
1199
1200        low and igh are integers determined by the balancing
1201          subroutine  balanc.  if  balanc  has not been used,
1202          set low=1, igh=n.
1203
1204        h contains the upper hessenberg matrix.
1205
1206        z contains the transformation matrix produced by  eltran
1207          after the reduction by  elmhes, or by  ortran  after the
1208          reduction by  orthes, if performed.  if the eigenvectors
1209          of the hessenberg matrix are desired, z must contain the
1210          identity matrix.
1211
1212        job  has the decimal decomposition xy;
1213            if x=0 hqror2 compute eigen-decomposition of h
1214            if x=1 hqror2 computes schur decomposition of h
1215            if x=2 eigenvalues are computed via schur decomposition
1216            if y=0 coordinate transformation is not accumulated
1217            if y=1 coordinate transformation is accumulated
1218
1219
1220     on output
1221
1222        h contains the Schur form
1223
1224        wr and wi contain the real and imaginary parts,
1225          respectively, of the eigenvalues.  the eigenvalues
1226          are unordered except that complex conjugate pairs
1227          of values appear consecutively with the eigenvalue
1228          having the positive imaginary part first.  if an
1229          error exit is made, the eigenvalues should be correct
1230          for indices ierr+1,...,n.
1231
1232        z contains the orthogonal transformation to the real schur
1233          form. If an error exit is made, z may be incorrect.
1234
1235        ierr is set to
1236          zero       for normal return,
1237          j          if the limit of 30*n iterations is exhausted
1238                     while the j-th eigenvalue is being sought.
1239
1240     calls cdiv for complex division.
1241
1242     questions and comments should be directed to burton s. garbow,
1243     mathematics and computer science div, argonne national laboratory
1244
1245     this version dated august 1983.
1246 */
1247 int dhqror2s(int _iLead, int _iSize, int _iLow, int _iHigh,
1248              double *_pdblHessUp, double *_pdblEigenReal, double *_pdblEigenImg,
1249              double *_pdblTrans, int _iMode)
1250 {
1251     //some loop ...
1252     int iIndex1         = 0;
1253     int iIndex2         = 0;
1254     int iIndex3         = 0;
1255     int iIndex4         = 0;
1256     int iIndex5         = 0;
1257     int iIndex6         = 0;
1258     int iIndex7         = 0;
1259     int iIndex8         = 0;
1260     int iIndex9         = 0;
1261     int iIndex10        = 0;
1262     int iIndex11        = 0;
1263     int iIndex12        = 0;
1264     int iIndex13        = 0;
1265     int iIndex14        = 0;
1266     int iIndex15        = 0;
1267     int iIndex16        = 0;
1268     int iIndex17        = 0;
1269     int iIndex18        = 0;
1270     int iIndex19        = 0;
1271     int iIndex20        = 0;
1272     int iIndex21        = 0;
1273     int iIndex22        = 0;
1274     int iIndex23        = 0;
1275     int iIndex24        = 0;
1276     int iIndex25        = 0;
1277     int iIndex26        = 0;
1278     int iIndex27        = 0;
1279     int iIndex28        = 0;
1280     int iIndex29        = 0;
1281     int iIndex30        = 0;
1282     int iIndex31        = 0;
1283     int iIndex32        = 0;
1284     int iIndex33        = 0;
1285
1286     int iLoop           = 0;
1287     int iPow10          = 0;
1288     int iModulo         = 0;
1289     int iOffset         = 0;
1290     int iMaxLoop        = 0;
1291
1292     int iIts            = 0;
1293     int iCoord1         = iOffset - 1;
1294     int iCoord2         = iCoord1 - 1;
1295     int iCoord3         = 0;
1296     int iLow            = _iLow - 1;
1297     int iHigh           = _iHigh - 1;
1298
1299     double dblNorm      = 0;
1300     double dblEps = nc_eps_machine();
1301
1302     double dblP         = 0;
1303     double dblQ         = 0;
1304     double dblR         = 0;
1305     double dblS         = 0;
1306     double dblT         = 0;
1307     double dblW         = 0;
1308     double dblX         = 0;
1309     double dblY         = 0;
1310     double dblZZ        = 0;
1311
1312     //I keep old source code, just for fun :
1313     //jx=job/10
1314     //jy=job-10*jx
1315     iPow10                      = _iMode / 10;
1316     iModulo                     = _iMode % 10;
1317
1318     //.......... store roots isolated by balanc
1319     //                          and compute matrix norm ..........
1320     for (iIndex1 = 0 ; iIndex1 < _iSize ; iIndex1++)
1321     {
1322         for (iIndex2 = iLoop ; iIndex2 < _iSize ; iIndex2++)
1323         {
1324             dblNorm += dabss(_pdblHessUp[iIndex1 + iIndex2 * _iSize]);
1325         }
1326
1327         iLoop = iIndex1;
1328         if (iPow10 == 1)
1329         {
1330             continue;
1331         }
1332
1333         if (iIndex1 >= iLow && iIndex1 <= iHigh)
1334         {
1335             continue;
1336         }
1337
1338         _pdblEigenReal[iIndex1] = _pdblHessUp[iIndex1 + iIndex1 * _iSize];
1339         _pdblEigenImg[iIndex1]  = 0;
1340     }
1341
1342     iOffset             = iHigh;
1343     dblT                = 0;
1344     iMaxLoop    = 30 * _iSize;
1345
1346     //.......... search for next eigenvalues ..........
1347 L60:
1348     if (iOffset <= iLow )
1349     {
1350         goto L340;
1351     }
1352
1353     iIts                = 0;
1354     iCoord1             = iOffset - 1;
1355     iCoord2             = iCoord1;
1356
1357 L70:
1358     //.......... look for single small sub-diagonal element
1359     //                  for l=en step -1 until low do -- ..........
1360     for (iIndex3 = iLow ; iIndex3 < (iOffset + 1); iIndex3++)
1361     {
1362         double dblTemp1 = 0, dblTemp2 = 0;
1363         iIndex4 = iOffset + iLow - iIndex3;
1364         if (iIndex4 == iLow)
1365         {
1366             break;
1367         }
1368
1369         dblS    = dabss(_pdblHessUp[(iIndex4 - 1) + (iIndex4 - 1) * _iSize]) +
1370                   dabss(_pdblHessUp[iIndex4 + iIndex4 * _iSize]);
1371
1372         if (dblS == 0)
1373         {
1374             dblS = dblNorm;
1375         }
1376
1377         dblTemp1        = dblS;
1378         dblTemp2        = dblTemp1 + dabss(_pdblHessUp[iIndex4 + (iIndex4 - 1) * _iSize]);
1379         if (dblTemp1 == dblTemp2)
1380         {
1381             break;
1382         }
1383     }
1384
1385     //.......... form shift ..........
1386     dblX        = _pdblHessUp[iOffset + iOffset * _iSize];
1387     if (iIndex4 == iOffset)
1388     {
1389         goto L270;
1390     }
1391
1392     dblY        = _pdblHessUp[iCoord1 + iCoord1 * _iSize];
1393     dblW        = _pdblHessUp[iOffset + iCoord1 * _iSize] * _pdblHessUp[iCoord1 + iOffset * _iSize];
1394     if (iIndex4 == iCoord1)
1395     {
1396         goto L280;
1397     }
1398
1399     if (iMaxLoop == 0)
1400     {
1401         goto L1000;
1402     }
1403
1404     if (iIts == 10 || iIts == 20)
1405     {
1406         //.......... form exceptional shift ..........
1407         dblT += dblX;
1408
1409         for (iIndex5 = iLow ; iIndex5 < iOffset ; iIndex5++)
1410         {
1411             _pdblHessUp[iIndex5 + iIndex5 * _iSize]     -= dblX;
1412         }
1413
1414         dblS    = dabss(_pdblHessUp[iOffset + iCoord1 * _iSize]) + dabss(_pdblHessUp[iCoord1 + iCoord2 * _iSize]);
1415         dblX    = 0.75 * dblS;
1416         dblY    = dblX;
1417         dblW    = -0.4375 * dblS * dblS;
1418     }
1419
1420
1421     iIts++;
1422     iMaxLoop--;
1423     //..........look for two consecutive small
1424     //                  sub-diagonal elements.
1425     //                  for m=en-2 step -1 until l do -- ..........
1426     for (iIndex6 = iIndex4; iIndex6 < iCoord2; iIndex6++)
1427     {
1428         double dblTemp1 = 0, dblTemp2 = 0;
1429
1430         iIndex7 = iCoord2 + (iIndex4 - 1) - iIndex6; // +/- 1
1431
1432         dblZZ   = _pdblHessUp[iIndex7 + iIndex7 * _iSize];
1433         dblR    = dblX - dblZZ;
1434         dblS    = dblY - dblZZ;
1435         dblP    = (dblR * dblS - dblW) / _pdblHessUp[(iIndex7 + 1) + iIndex7 * _iSize] + _pdblHessUp[iIndex7 + (iIndex7 + 1) * _iSize];
1436         dblQ    = _pdblHessUp[(iIndex7 + 1) + (iIndex7 + 1) * _iSize] - dblZZ - dblR - dblS;
1437         dblR    = _pdblHessUp[(iIndex7 + 2) + (iIndex7 + 1) * _iSize];
1438         dblS    = dabss(dblP) + dabss(dblQ) + dabss(dblR);
1439         dblP    /= dblS;
1440         dblQ    /= dblS;
1441         dblR    /= dblS;
1442
1443         if (iIndex7 == iIndex4)
1444         {
1445             break;
1446         }
1447
1448         dblTemp1 =      dabss(dblP) *
1449                     (dabss(_pdblHessUp[(iIndex7 - 1) + (iIndex7 - 1) * _iSize]) + dabss(dblZZ)) +
1450                     dabss(_pdblHessUp[(iIndex7 + 1) + (iIndex7 + 1) * _iSize]);
1451
1452         dblTemp2 =      dblTemp1 +
1453                     dabss(_pdblHessUp[iIndex7 + (iIndex7 - 1) * _iSize]) * (dabss(dblQ) + dabss(dblR));
1454         if (dblTemp1 == dblTemp2)
1455         {
1456             break;
1457         }
1458     }
1459     iCoord3     = iIndex7 + 2;
1460
1461     for (iIndex8 = iCoord3 ; iIndex8 <= iOffset ; iIndex8++)
1462     {
1463         _pdblHessUp[iIndex8 + (iIndex8 - 2) * _iSize]   = 0;
1464         if (iIndex8 != iCoord3)
1465         {
1466             _pdblHessUp[iIndex8 + (iIndex8 - 3) * _iSize]       = 0;
1467         }
1468     }
1469
1470     //..........        double qr step involving rows l to en and
1471     //                          columns m to en ..........
1472     for (iIndex9 = iIndex7 ; iIndex9 <= iCoord1 ; iIndex9++)
1473     {
1474         int iNotLast = iIndex9 != iCoord1;
1475         if (iIndex9 != iIndex7)
1476         {
1477             dblP                = _pdblHessUp[iIndex9 + (iIndex9 - 1) * _iSize];
1478             dblQ                = _pdblHessUp[(iIndex9 + 1) + (iIndex9 - 1) * _iSize];
1479             dblR                = 0;
1480             if (iNotLast)
1481             {
1482                 dblR    = _pdblHessUp[(iIndex9 + 2) + (iIndex9 - 1) * _iSize];
1483             }
1484             dblX                = dabss(dblP) + dabss(dblQ) + dabss(dblR);
1485             if (dblX == 0)
1486             {
1487                 continue;
1488             }
1489
1490             dblP        /= dblX;
1491             dblQ        /= dblX;
1492             dblR        /= dblX;
1493         }
1494
1495         dblS = dsigns(dsqrts(dblP * dblP + dblQ * dblQ + dblR * dblR), dblP);
1496
1497         if (iIndex9 != iIndex7)
1498         {
1499             _pdblHessUp[iIndex9 + (iIndex9 - 1) * _iSize] = -dblS * dblX;
1500         }
1501         else if (iIndex4 == iIndex7)
1502         {
1503             _pdblHessUp[iIndex9 + (iIndex9 - 1) * _iSize] *= -1;
1504         }
1505
1506         dblP    += dblS;
1507         dblX    = dblP / dblS;
1508         dblY    = dblQ / dblS;
1509         dblZZ   = dblR / dblS;
1510         dblQ    /= dblP;
1511         dblR    /= dblP;
1512
1513         if (iNotLast == 0)
1514         {
1515             //.......... row modification ..........
1516             for (iIndex10 = iIndex9 ; iIndex10 < _iSize ; iIndex10++)
1517             {
1518                 dblP    =       _pdblHessUp[iIndex9 + iIndex10 * _iSize] +
1519                             dblQ * _pdblHessUp[(iIndex9 + 1) + iIndex10 * _iSize];
1520
1521                 _pdblHessUp[iIndex9 + iIndex10 * _iSize]                -= dblP * dblX;
1522                 _pdblHessUp[(iIndex9 + 1) + iIndex10 * _iSize]  -= dblP * dblY;
1523             }
1524
1525             iIndex10 = Min(iOffset, iIndex9 + 3);
1526
1527             //.......... column modification ..........
1528             for (iIndex11 = 0 ; iIndex11 <= iIndex10 ; iIndex11++)
1529             {
1530                 dblP    =       dblX * _pdblHessUp[iIndex11 + iIndex9 * _iSize] +
1531                             dblY * _pdblHessUp[iIndex11 + (iIndex9 + 1) * _iSize];
1532
1533                 _pdblHessUp[iIndex11 + iIndex9 * _iSize]                -= dblP;
1534                 _pdblHessUp[iIndex11 + (iIndex9 + 1) * _iSize]  -= dblP * dblQ;
1535             }
1536
1537             if (iModulo == 1)
1538             {
1539                 //.......... accumulate transformations ..........
1540                 for (iIndex12 = iLow ; iIndex12 <= iHigh ; iIndex12++)
1541                 {
1542                     dblP        =       dblX * _pdblTrans[iIndex12 + iIndex9 * _iSize] +
1543                                 dblY * _pdblTrans[iIndex12 + (iIndex9 + 1) * _iSize];
1544
1545                     _pdblTrans[iIndex12 + iIndex9 * _iSize]                     -= dblP;
1546                     _pdblTrans[iIndex12 + (iIndex9 + 1) * _iSize]       -= dblP * dblQ;
1547                 }
1548             }
1549         }
1550         else
1551         {
1552             //.......... row modification ..........
1553             for (iIndex13 = iIndex9 ; iIndex13 < _iSize ; iIndex13++)
1554             {
1555                 dblP    =       _pdblHessUp[iIndex9 + iIndex13 * _iSize] +
1556                             dblQ * _pdblHessUp[(iIndex9 + 1) + iIndex13 * _iSize] +
1557                             dblR * _pdblHessUp[(iIndex9 + 2) + iIndex13 * _iSize];
1558
1559                 _pdblHessUp[iIndex9 + iIndex13 * _iSize]                -= dblP * dblX;
1560                 _pdblHessUp[(iIndex9 + 1) + iIndex13 * _iSize]  -= dblP * dblY;
1561                 _pdblHessUp[(iIndex9 + 2) + iIndex13 * _iSize]  -= dblP * dblZZ;
1562
1563             }
1564
1565             iIndex14 = Min(iOffset, iIndex9 + 3);
1566
1567             //.......... column modification ..........
1568             for (iIndex15 = 0 ; iIndex15 <= iIndex14 ; iIndex15++)
1569             {
1570                 dblP    =       dblX * _pdblHessUp[iIndex15 + iIndex9 * _iSize] +
1571                             dblY * _pdblHessUp[iIndex15 + (iIndex9 + 1) * _iSize] +
1572                             dblZZ * _pdblHessUp[iIndex15 + (iIndex9 + 2) * _iSize];
1573
1574                 _pdblHessUp[iIndex15 + iIndex9 * _iSize]                -= dblP;
1575                 _pdblHessUp[iIndex15 + (iIndex9 + 1) * _iSize]  -= dblP * dblQ;
1576                 _pdblHessUp[iIndex15 + (iIndex9 + 2) * _iSize]  -= dblP * dblR;
1577
1578             }
1579
1580             if (iModulo == 1)
1581             {
1582                 //.......... accumulate transformations ..........
1583                 for (iIndex16 = iLow ; iIndex16 <= iHigh ; iIndex16++)
1584                 {
1585                     dblP        =       dblX * _pdblTrans[iIndex16 + iIndex9 * _iSize] +
1586                                 dblY * _pdblTrans[iIndex16 + (iIndex9 + 1) * _iSize] +
1587                                 dblZZ * _pdblTrans[iIndex16 + (iIndex9 + 2) * _iSize];
1588
1589                     _pdblTrans[iIndex16 + iIndex9 * _iSize]                     -= dblP;
1590                     _pdblTrans[iIndex16 + (iIndex9 + 1) * _iSize]       -= dblP * dblQ;
1591                     _pdblTrans[iIndex16 + (iIndex9 + 2) * _iSize]       -= dblP * dblR;
1592                 }
1593             }
1594         }
1595     }
1596
1597     goto L70;
1598
1599     //.......... one root found ..........
1600 L270:
1601     _pdblHessUp[iOffset + iOffset * _iSize] = dblX + dblT;
1602     //ADDED TO MARK BLOCK SEPARATION BY HARD ZEROS
1603     if (iOffset + 1 < _iSize)
1604     {
1605         _pdblHessUp[(iOffset + 1) + (iOffset + 1) * _iSize] = 0;
1606     }
1607
1608     if (iPow10 != 1)
1609     {
1610         _pdblEigenReal[iOffset] = _pdblHessUp[iOffset + iOffset * _iSize];
1611         _pdblEigenImg[iOffset]  = 0;
1612     }
1613     iOffset = iCoord1;
1614     goto L60;
1615
1616     //.......... two roots found ..........
1617 L280:
1618     dblP                                                                        = (dblY - dblX) / 2.0;
1619     dblQ                                                                        = dblP * dblP + dblW;
1620     dblZZ                                                                       = dsqrts(dabss(dblQ));
1621     _pdblHessUp[iOffset + iOffset * _iSize] = dblX + dblT;
1622     dblX                                                                        = _pdblHessUp[iOffset + iOffset * _iSize];
1623     _pdblHessUp[iCoord1 + iCoord1 * _iSize]     = dblY + dblT;
1624
1625
1626     if (dblQ < 0)
1627     {
1628         goto L320;
1629     }
1630
1631     //.......... real pair ..........
1632     dblZZ       = dblP + dsigns(dblZZ, dblP);
1633     if (iPow10 != 1)
1634     {
1635         _pdblEigenReal[iCoord1]         = dblX + dblZZ;
1636         _pdblEigenReal[iOffset]         = _pdblEigenReal[iCoord1];
1637
1638         if (dblZZ != 0)
1639         {
1640             _pdblEigenReal[iOffset]     = dblX - dblW / dblZZ;
1641         }
1642
1643         _pdblEigenImg[iCoord1]          = 0;
1644         _pdblEigenImg[iOffset]          = 0;
1645     }
1646
1647     dblX        = _pdblHessUp[iOffset + iCoord1 * _iSize];
1648     dblS        = dabss(dblX) + dabss(dblZZ);
1649     dblP        = dblX / dblS;
1650     dblQ        = dblZZ / dblS;
1651     dblR        = dsqrts(dblP * dblP + dblQ * dblQ);
1652     dblP        = dblP / dblR;
1653     dblQ        = dblQ / dblR;
1654
1655     //.......... row modification ..........
1656     for (iIndex17 = iCoord1 ; iIndex17 < _iSize ; iIndex17++)
1657     {
1658         dblZZ   = _pdblHessUp[iCoord1 + iIndex17 * _iSize];
1659
1660         _pdblHessUp[iCoord1 + iIndex17 * _iSize]        = dblQ * dblZZ + dblP * _pdblHessUp[iOffset + iIndex17 * _iSize];
1661         _pdblHessUp[iOffset + iIndex17 * _iSize]        = dblQ * _pdblHessUp[iOffset + iIndex17 * _iSize] - dblP * dblZZ;
1662     }
1663
1664     //.......... column modification ..........
1665     for (iIndex18 = 0 ; iIndex18 <= iOffset ; iIndex18++)
1666     {
1667         dblZZ   = _pdblHessUp[iIndex18 + iCoord1 * _iSize];
1668
1669         _pdblHessUp[iIndex18 + iCoord1 * _iSize]        = dblQ * dblZZ + dblP * _pdblHessUp[iIndex18 + iOffset * _iSize];
1670         _pdblHessUp[iIndex18 + iOffset * _iSize]        = dblQ * _pdblHessUp[iIndex18 + iOffset * _iSize] - dblP * dblZZ;
1671     }
1672
1673     //.......... accumulate transformations ..........
1674     if (iModulo == 1)
1675     {
1676         for (iIndex19 = iLow ; iIndex19 <= iHigh ; iIndex19++)
1677         {
1678             dblZZ       = _pdblTrans[iIndex19 + iCoord1 * _iSize];
1679
1680             _pdblTrans[iIndex19 + iCoord1 * _iSize]     = dblQ * dblZZ + dblP * _pdblTrans[iIndex19 + iOffset * _iSize];
1681             _pdblTrans[iIndex19 + iOffset * _iSize]     = dblQ * _pdblTrans[iIndex19 + iOffset * _iSize] - dblP * dblZZ;
1682         }
1683     }
1684
1685     //ADDED TO MARK BLOCK SEPARATION BY HARD ZEROS
1686     _pdblHessUp[iOffset + iCoord1 * _iSize]                     = 0;
1687     if (iOffset + 1 != _iSize)
1688     {
1689         _pdblHessUp[iOffset + 1 + iOffset * _iSize] = 0;
1690     }
1691
1692     goto L330;
1693
1694     //.......... complex pair ..........
1695 L320:
1696     if (iPow10 != 1)
1697     {
1698         _pdblEigenReal[iCoord1] = dblX + dblP;
1699         _pdblEigenReal[iOffset] = dblX + dblP;
1700         _pdblEigenImg[iCoord1]  = dblZZ;
1701         _pdblEigenImg[iOffset]  = -dblZZ;
1702     }
1703
1704     //ADDED TO MARK BLOCK SEPARATION BY HARD ZEROS
1705     if (iOffset + 1 <= _iSize)
1706     {
1707         _pdblHessUp[iOffset + 1 + iOffset * _iSize] = 0;
1708     }
1709
1710 L330:
1711     iOffset     = iCoord2;
1712     goto L60;
1713
1714 L340:
1715     if (iPow10 != 0 || dblNorm == 0)
1716     {
1717         goto L1001;
1718     }
1719
1720     //:::::::::: for en=n step -1 until 1 do -- ::::::::::
1721     for (iIndex21 = 0 ; iIndex21 < _iSize ; iIndex21++)
1722     {
1723         double dblCRES  = 0;
1724
1725         iOffset = _iSize + 1 - iIndex21;
1726         dblP    = _pdblEigenReal[iOffset];
1727         dblQ    = _pdblEigenImg[iOffset];
1728         iCoord1 = iOffset - 1;
1729         dblQ    = dblQ + 1;
1730         dblCRES = dblQ - 1;
1731         if (dblCRES < 0)
1732         {
1733             goto L710;
1734         }
1735         else if (dblCRES > 0)
1736         {
1737             continue;
1738         }
1739
1740         //:::::::::: real vector ::::::::::
1741         iIndex22        = iOffset;
1742         _pdblHessUp[iOffset + iOffset * _iSize] = 1;
1743         if (iCoord1 == 0)
1744         {
1745             continue;
1746         }
1747
1748         //:::::::::: for i=en-1 step -1 until 1 do -- ::::::::::
1749         for (iIndex23 = 0 ; iIndex23 < iCoord1 ; iIndex23++)
1750         {
1751             iIndex24    = iOffset - iIndex23;
1752             dblW                = _pdblHessUp[iIndex24 + iIndex24 * _iSize] - dblP;
1753             dblR                = _pdblHessUp[iIndex24 + iOffset * _iSize];
1754             if (iIndex22 <= iCoord1)
1755                 for (iIndex25 = iIndex22 ; iIndex25 < iCoord1 ; iIndex25++)
1756                     dblR += _pdblHessUp[iIndex24 + iIndex25 * _iSize] *
1757                             _pdblHessUp[iIndex25 + iOffset * _iSize];
1758
1759             if (_pdblEigenImg[iIndex24] >= 0)
1760             {
1761                 goto L630;
1762             }
1763             dblZZ       = dblW;
1764             dblS        = dblR;
1765             continue;
1766
1767 L630:
1768             iIndex22    = iIndex24;
1769             if (_pdblEigenImg[iIndex24] != 0)
1770             {
1771                 goto L640;
1772             }
1773
1774             dblT                = dblW;
1775             if (dblW == 0)
1776             {
1777                 dblT    = dblEps * dblNorm;
1778             }
1779
1780             _pdblHessUp[iIndex24 + iOffset * _iSize] = -dblR / dblT;
1781             continue;
1782
1783 L640:
1784             //:::::::::: solve real equations ::::::::::
1785             dblX        = _pdblHessUp[iIndex24 + (iIndex24 + 1) * _iSize];
1786             dblY        = _pdblHessUp[(iIndex24 + 1) + iIndex24 * _iSize];
1787             dblQ        =       (_pdblEigenReal[iIndex24] - dblP) *
1788                         (_pdblEigenReal[iIndex24] - dblP) +
1789                         _pdblEigenImg[iIndex24] * _pdblEigenImg[iIndex24];
1790             dblT        = (dblX * dblS - dblZZ * dblR) / dblQ;
1791             _pdblHessUp[iIndex24 + iOffset * _iSize]    = dblT;
1792
1793             if (dabss(dblX) > dabss(dblZZ))
1794             {
1795                 _pdblHessUp[iIndex24 + 1 + iOffset * _iSize] = (-dblR - dblX * dblW) / dblX;
1796             }
1797             else
1798             {
1799                 _pdblHessUp[iIndex24 + 1 + iOffset * _iSize] = (-dblR - dblX * dblW) / dblX;
1800             }
1801         }
1802         //:::::::::: end real vector ::::::::::
1803         continue;
1804         //:::::::::: complex vector ::::::::::
1805 L710:
1806         iIndex22        = iCoord1;
1807         //::::::::::    last vector component chosen imaginary so that
1808         //                              eigenvector matrix is triangular ::::::::::
1809         if (dabss(_pdblHessUp[iOffset + iCoord1 * _iSize]) > dabss(_pdblHessUp[iCoord1 + iOffset * _iSize]))
1810         {
1811             _pdblHessUp[iCoord1 + iCoord1 * _iSize]     =       dblQ / _pdblHessUp[iOffset + iCoord1 * _iSize];
1812             _pdblHessUp[iCoord1 + iOffset * _iSize]     =       -(_pdblHessUp[iOffset + iOffset * _iSize] - dblP) /
1813                     _pdblHessUp[iOffset + iCoord1 * _iSize];
1814         }
1815         else
1816         {
1817             double dblReal      = _pdblHessUp[iCoord1 + iCoord1 * _iSize] - dblP; //z3r
1818             double dblImg       = 0; //z3i
1819             double dblCoef      = dblReal * dblReal + dblQ * dblQ; //z3
1820
1821             _pdblHessUp[iCoord1 + iCoord1 * _iSize]     =
1822                 -_pdblHessUp[iCoord1 + iOffset * _iSize] * dblQ / dblCoef;
1823             _pdblHessUp[iCoord1 + iOffset * _iSize]     =
1824                 -_pdblHessUp[iCoord1 + iOffset * _iSize] * dblReal / dblCoef;
1825             _pdblHessUp[iOffset + iCoord1 * _iSize]     = 0;
1826             _pdblHessUp[iOffset + iOffset * _iSize]     = 1;
1827             iCoord2     = iCoord1 - 1;
1828
1829             if (iCoord2 == 0)
1830             {
1831                 continue;
1832             }
1833
1834             //:::::::::: for i=en-2 step -1 until 1 do -- ::::::::::
1835             for (iIndex25 = 0 ; iIndex25 < iCoord2 ; iIndex25++)
1836             {
1837                 double dblTemp11        = 0; //r
1838                 double dblTemp12        = 0; //ra
1839                 double dblTemp21        = 0; //s
1840                 double dblTemp22        = 0; //sa
1841
1842                 iIndex26                = iCoord1 - iIndex25;
1843                 dblW                    = _pdblHessUp[iIndex24 + iIndex24 * _iSize] - dblP;
1844                 dblTemp12               = 0;
1845                 dblTemp22               = _pdblHessUp[iIndex24 + iOffset * _iSize];
1846
1847                 for (iIndex27 = iIndex22 ; iIndex27 < iCoord1 ; iIndex27++)
1848                 {
1849                     dblTemp12   +=      _pdblHessUp[iIndex24 + iIndex27 * _iSize] *
1850                                     _pdblHessUp[iIndex27 + iCoord1 * _iSize];
1851                     dblTemp22   +=      _pdblHessUp[iIndex24 + iIndex27 * _iSize] *
1852                                     _pdblHessUp[iIndex27 + iOffset * _iSize];
1853                 }
1854
1855                 if (_pdblEigenImg[iIndex24] < 0)
1856                 {
1857                     dblZZ               = dblW;
1858                     dblTemp11   = dblTemp12;
1859                     dblTemp21   = dblTemp22;
1860                 }
1861                 else
1862                 {
1863
1864                     if (_pdblEigenImg[iIndex24] == 0)
1865                     {
1866                         dblCoef =       dblW * dblW + dblQ * dblQ;
1867                         dblReal = -dblTemp12 * dblW - dblTemp22 * dblQ;
1868                         dblImg  = dblTemp12 * dblQ - dblTemp22 * dblW;
1869                         _pdblHessUp[iIndex24 + iCoord1 * _iSize] = dblReal / dblCoef;
1870                         _pdblHessUp[iIndex24 + iOffset * _iSize] = dblImg / dblCoef;
1871                     }
1872                     else
1873                     {
1874                         double dblVectReal      = 0;
1875                         double dblVectImg       = 0;
1876
1877                         dblX                            =       _pdblHessUp[iIndex24 + (iIndex24 + 1) * _iSize];
1878                         dblY                            =       _pdblHessUp[(iIndex24 + 1) + iIndex24 * _iSize];
1879                         dblVectReal                     =       (_pdblEigenReal[iIndex24] - dblP) * (_pdblEigenReal[iIndex24] - dblP) +
1880                                                 _pdblEigenImg[iIndex24] * _pdblEigenImg[iIndex24] - dblQ * dblQ;
1881                         dblVectImg                      =       (_pdblEigenReal[iIndex24] - dblP) * 2 * dblQ;
1882
1883                         if (dblVectReal == 0 && dblVectImg == 0)
1884                             dblVectReal         =       dblEps * dblNorm *
1885                                                 (dabss(dblW) + dabss(dblQ) + dabss(dblX) + dabss(dblY) + dabss(dblZZ));
1886
1887                         dblReal                         = dblX * dblR - dblZZ * dblTemp12 + dblQ * dblTemp22;
1888                         dblImg                          = dblX * dblS - dblZZ * dblTemp22 - dblQ * dblTemp12;
1889                         dblCoef                         = dblVectReal * dblVectReal + dblVectImg * dblVectImg;
1890
1891                         _pdblHessUp[iIndex24 + iCoord1 * _iSize]        = (dblReal * dblVectReal + dblImg * dblVectImg) / dblCoef;
1892                         _pdblHessUp[iIndex24 + iOffset * _iSize]        = (-dblReal * dblVectImg + dblImg * dblVectReal) / dblCoef;
1893
1894                         if (dabss(dblX) > dabss(dblZZ) + dabss(dblQ))
1895                         {
1896                             _pdblHessUp[iIndex24 + 1 + iCoord1 * _iSize] =
1897                                 (-dblTemp12 - dblW * _pdblHessUp[iIndex24 + iCoord1 * _iSize] + dblQ * _pdblHessUp[iIndex24 + iOffset * _iSize]) / dblX;
1898                             _pdblHessUp[iIndex24 + 1 + iOffset * _iSize] =
1899                                 (-dblTemp22 - dblW * _pdblHessUp[iIndex24 + iOffset * _iSize] + dblQ * _pdblHessUp[iIndex24 + iCoord1 * _iSize]) / dblX;
1900                         }
1901                         else
1902                         {
1903                             dblReal     = -dblR - dblY * _pdblHessUp[iIndex24 + iCoord1 * _iSize];
1904                             dblImg      = -dblS - dblY * _pdblHessUp[iIndex24 + iOffset * _iSize];
1905                             dblCoef     = dblZZ * dblZZ + dblQ * dblQ;
1906                             _pdblHessUp[iIndex24 + 1 + iCoord1 * _iSize] = (dblReal * dblZZ + dblImg * dblQ) / dblCoef;
1907                             _pdblHessUp[iIndex24 + 1 + iOffset * _iSize] = (-dblReal * dblQ + dblImg * dblZZ) / dblCoef;
1908                         }
1909                     }
1910                 }
1911             }
1912             //:::::::::: end complex vector ::::::::::
1913         }
1914         //:::::::::: end back substitution.
1915     }
1916
1917     if (iModulo == 0)
1918     {
1919         return 0;
1920     }
1921
1922     //vectors of isolated roots ::::::::::
1923     for (iIndex27 = 0 ; iIndex27 < _iSize ; iIndex27++)
1924     {
1925         if (iIndex27 >= iLow && iIndex27 <= iHigh)
1926         {
1927             continue;
1928         }
1929
1930         for (iIndex28 = iIndex27 ; iIndex28 < _iSize ; iIndex28++)
1931         {
1932             _pdblTrans[iIndex27 + iIndex28 * _iSize] = _pdblHessUp[iIndex27 + iIndex28 * _iSize];
1933         }
1934     }
1935
1936     //:::::::::: multiply by transformation matrix to give
1937     //                  vectors of original full matrix.
1938     //                  for j=n step -1 until low do -- ::::::::::
1939     for (iIndex29 = iLow ; iIndex29 < _iSize ; iIndex29++)
1940     {
1941         for (iIndex32 = iLow ; iIndex32 < iHigh ; iIndex32++)
1942         {
1943             dblZZ = 0;
1944             for (iIndex33 = iLow ; iIndex33 < iIndex31 ; iIndex33++)
1945             {
1946                 dblZZ += _pdblTrans[iIndex32 + iIndex33 * _iSize] * _pdblHessUp[iIndex33 + iIndex30 * _iSize];
1947                 _pdblTrans[iIndex32 + iIndex30 * _iSize] = dblZZ;
1948             }
1949         }
1950     }
1951     //..........        set error -- all eigenvalues have not
1952     //                          converged after 30*n iterations ..........
1953 L1001:
1954     return iOffset;
1955 L1000:
1956     return 0;
1957 }
1958 /*--------------------------------------------------------------------------*/
1959
1960 int dpades(double *_pdblVal, int _iLeadDimIn, int _iSize, double *_pdblExp, int _iLeadDimOut, double *_pdblAlpha, double *_pdblWS, int *_piWS)
1961 {
1962     int iErr            = 0;
1963
1964     double dblZero      = 0;
1965     double dblOne       = 1;
1966     double dblTwo       = 2;
1967     double dblEfact     = 0;
1968     double dblNorm      = 0;
1969     double dblRcon      = 0;
1970
1971     int iZero   = 0;
1972     int iMax    = 10;
1973     int iIndex1 = 0;
1974     int iIndex2 = 0;
1975     int iLoop1  = 0;
1976     int iLoop2  = 0;
1977     int iSize2  = _iSize * _iSize;
1978
1979     if (sdblExpmN < 0)
1980     {
1981         //compute de pade's aprroximants type which is necessary to obtain
1982         //machine precision
1983         C2F(coef)(&iErr);
1984         if (iErr != 0)
1985         {
1986             return iErr;
1987         }
1988     }
1989
1990     iIndex1             = 0;
1991     dblEfact    = dblOne;
1992
1993     if (*_pdblAlpha < 1)
1994     {
1995         goto L90;
1996     }
1997
1998     for (iLoop1 = 0 ; iLoop1 < iMax ; iLoop1++)
1999     {
2000         iIndex1++;
2001         dblEfact *= dblTwo;
2002
2003         if (*_pdblAlpha <= dblEfact)
2004         {
2005             goto L60;
2006         }
2007     }
2008
2009 L30:
2010     iIndex1++;
2011     dblEfact *= dblTwo;
2012     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2013         for (iLoop2 = 0 ; iLoop2 < _iSize ; iLoop2++)
2014         {
2015             _pdblVal[iLoop1 + iLoop2 * _iSize] /= dblTwo;
2016         }
2017
2018     dblNorm /= dblTwo;
2019     goto L115;
2020
2021     //we find a matrix a'=a*2-m whith a spectral radius smaller than one.
2022 L60:
2023     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2024         for (iLoop2 = 0 ; iLoop2 < _iSize ; iLoop2++)
2025         {
2026             _pdblVal[iLoop1 + iLoop2 * _iSize] /= dblEfact;
2027         }
2028
2029 L90:
2030     C2F(cerr)(_pdblVal, _pdblWS, &_iLeadDimIn, &_iSize, &sdblExpmN, &iIndex1, &iMax);
2031
2032     dblNorm     = dblZero;
2033     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2034     {
2035         *_pdblAlpha = dblZero;
2036         for (iLoop2 = 0 ; iLoop2 < _iSize ; iLoop2++)
2037         {
2038             *_pdblAlpha += dabss(_pdblVal[iLoop1 + iLoop2 * _iSize]);
2039         }
2040         if (*_pdblAlpha > dblNorm)
2041         {
2042             dblNorm = *_pdblAlpha;
2043         }
2044     }
2045
2046     //compute the inverse of the denominator of dpade's approximants.
2047 L115:
2048     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2049         for (iLoop2 = 0 ; iLoop2 < _iSize ; iLoop2++)
2050         {
2051             _pdblExp[iLoop1 + iLoop2 * _iSize] = -_pdblVal[iLoop1 + iLoop2 * _iSize];
2052         }
2053
2054     C2F(dclmat)(&_iLeadDimOut, &_iSize, _pdblExp, _pdblWS, &_iSize, &_pdblWS[iSize2 + 1], spdblExpmC, &sdblExpmN);
2055
2056     //compute de l-u decomposition of n (-a) and the condition numbwk(n2+1)
2057     //                                                                  pp
2058     C2F(dgeco)(_pdblWS, &_iSize, &_iSize, _piWS, &dblRcon, &_pdblWS[iSize2 + 1]);
2059
2060     dblRcon = dblRcon * dblRcon * dblRcon * dblRcon;
2061
2062     if (dblRcon + dblOne <= dblOne && dblNorm > dblOne && iIndex1 < iMax)
2063     {
2064         goto L30;
2065     }
2066
2067     //compute the numerator of dpade's approximants.
2068     C2F(dclmat)(&_iLeadDimIn, &_iSize, _pdblVal, _pdblExp, &_pdblWS[iSize2 + 1], spdblExpmC, &sdblExpmN);
2069
2070     //compute the dpade's approximants by
2071
2072     //  n (-a) x=n (a)
2073     //  pp      pp
2074     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2075     {
2076         C2F(dgesl)(_pdblWS, &_iSize, &_iSize, _piWS, &_pdblExp[iLoop1 * _iSize], &iZero);
2077     }
2078
2079     if (iIndex1 != 0)
2080     {
2081         //remove the effects of normalization.
2082         for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2083         {
2084             C2F(dmmul)(_pdblExp, &_iLeadDimOut, _pdblExp, &_iLeadDimOut, _pdblWS, &_iSize, &_iSize, &_iSize, &_iSize);
2085             C2F(dmcopy)(_pdblWS, &_iSize, _pdblExp, &_iLeadDimOut, &_iSize, &_iSize);
2086         }
2087     }
2088
2089     return 0;
2090 }
2091
2092 /*
2093      PURPOSE
2094         computes the matrix product C = A * B
2095             C   =   A   *   B
2096           (l,n)   (l,m) * (m,n)
2097
2098      PARAMETERS
2099         input
2100         -----
2101         A : (double) array (l, m) with leading dim na
2102
2103         B : (double) array (m, n) with leading dim nb
2104
2105         na, nb, nc, l, m, n : integers
2106
2107         output
2108         ------
2109         C : (double) array (l, n) with leading dim nc
2110
2111      NOTE
2112         (original version substituted by a call to the blas dgemm)
2113 */
2114 void ddmmuls(double *_pdblA, int _iLeadDimA,
2115              double *_pdblB, int _iLeadDimB,
2116              double *_pdblOut, int _iLeadDimOut,
2117              int _iRowsA, int _iColsA, int _iColsB)
2118 {
2119     char cN                     = 'n';
2120     double dblOne       = 1;
2121     double dblZero      = 0;
2122
2123     C2F(dgemm)(&cN, &cN, &_iRowsA, &_iColsB, &_iColsA, &dblOne, _pdblA, &_iLeadDimA, _pdblB, &_iLeadDimB, &dblZero, _pdblOut, &_iLeadDimOut);
2124 }
2125
2126
2127 /*
2128         purpose
2129          given  upper hessenberg matrix a
2130          with consecutive ls1xls1 and ls2xls2 diagonal blocks (ls1,ls2.le.2)
2131          starting at row/column l, exch produces equivalence transforma-
2132          tion zt that exchange the blocks along with their
2133          eigenvalues.
2134
2135         syntax
2136
2137                 subroutine exch(nmax,n,a,z,l,ls1,ls2)
2138                 integer nmax,n,l,ls1,ls2
2139                 double precision a(nmax,n),z(nmax,n)
2140
2141                 nmax     the first dimension of a, b and z
2142                 n        the order of a, and z
2143            *a        the matrix whose blocks are to be interchanged
2144            *z        upon return this array is multiplied by the column
2145                                  transformation zt.
2146                 l        the position of the blocks
2147                 ls1      the size of the first block
2148                 ls2      the size of the second block
2149
2150         auxiliary routines
2151                 drot (blas)
2152                 giv
2153                 max abs (fortran)
2154         originator
2155                 Delebecque f. and Steer s. INRIA adapted from exchqz (VanDooren)
2156 */
2157
2158 int dexchs(int _iMax, int _iLeadDim, double *_pdblIn, double *_pdblOut,
2159            int _iPos, int _iSize1, int _iSize2)
2160 {
2161     int iLoop1          = 0; //i
2162     int iLoop2          = 0; //j
2163     int iIndex1         = 0; //l1
2164     int iIndex2         = 0; //l2
2165     int iIndex3         = 0; //l3
2166     int iIndex4         = 0; //li
2167     int iIndex5         = 0; //lj
2168     int iIndex6         = 0; //ll
2169
2170     //For Fortran compatibility
2171     int iZero           = 0;
2172     int iOne            = 1;
2173     int iTwo            = 2;
2174     int iThree          = 3;
2175     int iFour           = 4;
2176
2177
2178     double pdblTemp[3][4] = {0};
2179     double dblD         = 0;
2180     double dblE         = 0;
2181     double dblF         = 0;
2182     double dblG         = 0;
2183     double dblSa        = 0;
2184     double dblSb        = 0;
2185
2186     double dblTemp1     = 0;
2187     double dblTemp2     = 0;
2188     double dblTemp3     = 0;
2189
2190     iIndex1 = _iPos + 1;
2191     iIndex6     = _iSize1 + _iSize2;
2192
2193     if (iIndex6 <= 2)
2194     {
2195
2196         //interchange 1x1 and 1x1 blocks via an equivalence
2197         //transformation       a:=z'*a*z ,
2198         //where z is givens rotation
2199
2200         dblF    = Max(dabss(_pdblIn[(iIndex1 - 1) + (iIndex1 - 1) * _iLeadDim]), 1);
2201         dblSa   = _pdblIn[(iIndex1 - 1) + (iIndex1 - 1) * _iLeadDim] / dblF;
2202         dblSb   = 1.0 / dblF;
2203         dblF    = dblSa - dblSb * _pdblIn[_iPos + _iPos * _iLeadDim];
2204
2205         //construct the column transformation z
2206         dblG    = -dblSa * _pdblIn[_iPos + (iIndex1 - 1) * _iLeadDim];
2207         dgivs(dblF, dblG, &dblD, &dblE);
2208
2209         dblTemp1 = -dblD;
2210         C2F(drot)(&iIndex1, &_pdblIn[_iPos * _iLeadDim], &iOne, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2211         C2F(drot)(&_iLeadDim, &_pdblOut[_iPos * _iLeadDim], &iOne, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2212
2213         //construct the row transformation q
2214         dblTemp2 = _iLeadDim - _iPos + 1;
2215         C2F(drot)(&dblTemp2, &_pdblIn[_iPos + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &iOne, &dblE, &dblTemp1);
2216         _pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim]      = 0;
2217         return 0;
2218     }
2219
2220     //interchange 1x1 and 2x2 blocks via an equivalence
2221     //transformation  a:=z2'*z1'*a*z1*z2 ,
2222     //where each zi is a givens rotation
2223     iIndex2     = _iPos + 2;
2224
2225     if (_iSize1 != 2)
2226     {
2227         dblG    = Max(dabss(_pdblIn[_iPos + _iPos * _iLeadDim]), 1);
2228         //evaluate the pencil at the eigenvalue corresponding
2229         //to the 1x1 block
2230         dblSa   = _pdblIn[_iPos + _iPos * _iLeadDim] / dblG;
2231         dblSb   = 1 / dblG;
2232         for (iLoop1 = 0 ; iLoop1 < 2 ; iLoop1++)
2233         {
2234             iIndex5 = _iPos + iLoop1;
2235             for (iLoop2 = 0 ; iLoop2 < 3 ; iLoop2++)
2236             {
2237                 iIndex4                                         = _iPos + iLoop2 - 1;
2238                 pdblTemp[iLoop2][iLoop1]        = -dblSb * _pdblIn[iIndex4 + iIndex5 * _iLeadDim];
2239                 if (iLoop1 == iLoop2)
2240                 {
2241                     pdblTemp[iLoop2][iLoop1] += dblSa;
2242                 }
2243             }//for
2244         }//for
2245
2246         dgivs(pdblTemp[2][0], pdblTemp[2][1], &dblD, &dblE);
2247         dblTemp1 = -dblD;
2248         C2F(drot)(&iThree, &pdblTemp[0][0], &iOne, &pdblTemp[0][1], &iOne, &dblE, &dblTemp1);
2249         //perform the row transformation z1'
2250         dgivs(pdblTemp[0][0], pdblTemp[1][0], &dblD, &dblE);
2251         pdblTemp[1][1] = -pdblTemp[0][1] * dblE + pdblTemp[1][1] * dblD;
2252         dblTemp2 = _iLeadDim - _iPos + 1;
2253         C2F(drot)(&dblTemp2, &_pdblIn[_iPos + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &_iMax, &dblD, &dblE);
2254         //perform the column transformation z1
2255         C2F(drot)(&iIndex2, &_pdblIn[_iPos * _iLeadDim], &iOne, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2256         C2F(drot)(&_iLeadDim, &_pdblOut[_iPos * _iLeadDim], &iOne, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2257         //perform the row transformation z2'
2258         dgivs(pdblTemp[1][1], pdblTemp[2][1], &dblD, &dblE);
2259         C2F(drot)(&dblTemp2, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim], &_iMax, &dblD, &dblE);
2260         //perform the column transformation z2
2261         C2F(drot)(&iIndex2, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &_pdblIn[(iIndex2 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2262         C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &_pdblOut[(iIndex2 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2263         //put the neglectable elements equal to zero
2264         _pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim] = 0;
2265         _pdblIn[(iIndex2 - 1) + (iIndex1 - 1) * _iLeadDim] = 0;
2266         return 0;
2267     }
2268
2269     //interchange 2x2 and 1x1 blocks via an equivalence
2270     //transformation  a:=z2'*z1'*a*z1*z2 ,
2271     //where each zi is a givens rotation
2272     if (_iSize2 != 2)
2273     {
2274         dblG    = Max(dabss(_pdblIn[(iIndex2 - 1) + (iIndex2 - 1) * _iLeadDim]), 1);
2275         //evaluate the pencil at the eigenvalue corresponding
2276         //to the 1x1 block
2277         dblSa   = _pdblIn[(iIndex2 - 1) + (iIndex2 - 1) * _iLeadDim] / dblG;
2278         dblSb   = 1 / dblG;
2279         for (iLoop1 = 0 ; iLoop1 < 2 ; iLoop1++)
2280         {
2281             iIndex5 = _iPos + iLoop1 - 1;
2282             for (iLoop2 = 0 ; iLoop2 < 3 ; iLoop2++)
2283             {
2284                 iIndex4                                         = _iPos + iLoop2 - 1;
2285                 pdblTemp[iLoop1][iLoop2]        = -dblSb * _pdblIn[iIndex5 + iIndex4 * _iLeadDim];
2286                 if (iLoop1 == iLoop2)
2287                 {
2288                     pdblTemp[iLoop1][iLoop2] += dblSa;
2289                 }
2290             }//for
2291         }//for
2292
2293         dgivs(pdblTemp[0][0], pdblTemp[1][0], &dblD, &dblE);
2294         C2F(drot)(&iThree, &pdblTemp[0][0], &iThree, &pdblTemp[1][0], &iThree, &dblD, &dblE);
2295         //perform the row transformation z1
2296         dgivs(pdblTemp[1][1], pdblTemp[1][2], &dblD, &dblE);
2297         pdblTemp[0][1] = pdblTemp[0][1] * dblE - pdblTemp[0][2] * dblD;
2298         dblTemp1 = -dblD;
2299         C2F(drot)(&iIndex2, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &_pdblIn[(iIndex2 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2300         C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &_pdblOut[(iIndex2 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2301         //perform the column transformation z1'
2302         dblTemp2 = _iLeadDim - _iPos + 1;
2303         C2F(drot)(&dblTemp2, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim], &_iMax, &dblE, &dblTemp1);
2304         //perform the row transformation z2
2305         dgivs(pdblTemp[0][0], pdblTemp[0][1], &dblD, &dblE);
2306         C2F(drot)(&dblTemp2, &_pdblIn[_iPos * _iLeadDim], &iOne, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2307         C2F(drot)(&_iLeadDim, &_pdblOut[_iPos * _iLeadDim], &iOne, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &dblE, &dblTemp1);
2308         //perform the column transformation z2'
2309         C2F(drot)(&dblTemp2, &_pdblIn[_iPos + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &iOne, &dblE, &dblTemp1);
2310         //put the neglectable elements equal to zero
2311         _pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim] = 0;
2312         _pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim] = 0;
2313         return 0;
2314     }
2315
2316     iIndex3     = _iPos + 3;
2317     dblD        =       _pdblIn[(iIndex2 - 1) + (iIndex2 - 1) * _iLeadDim] *
2318                 _pdblIn[(iIndex3 - 1) + (iIndex3 - 1) * _iLeadDim] -
2319                 _pdblIn[(iIndex2 - 1) + (iIndex3 - 1) * _iLeadDim] *
2320                 _pdblIn[(iIndex3 - 1) + (iIndex2 - 1) * _iLeadDim];
2321
2322     dblE        =       _pdblIn[(iIndex2 - 1) + (iIndex2 - 1) * _iLeadDim] +
2323                 _pdblIn[(iIndex3 - 1) + (iIndex3 - 1) * _iLeadDim];
2324
2325     ddmmuls(    &_pdblIn[_iPos + _iPos * _iLeadDim], _iMax,
2326                 &_pdblIn[_iPos + _iPos * _iLeadDim], _iMax,
2327                 pdblTemp[0], 3, 2, 4, 4);
2328
2329     for (iLoop1 = 0 ; iLoop1 < 2 ; iLoop1++)
2330     {
2331         pdblTemp[iLoop1][iLoop1] += dblD;
2332         for (iLoop2 = 0 ; iLoop2 < 4 ; iLoop2++)
2333         {
2334             pdblTemp[iLoop1][iLoop2] = pdblTemp[iLoop1][iLoop2] - dblE *
2335                                        _pdblIn[(_iPos - 1 + iLoop1) + (_iPos - 1 + iLoop2) * _iLeadDim];
2336         }
2337     }
2338
2339     //g0
2340     dgivs(pdblTemp[0][0], pdblTemp[1][0], &dblD, &dblE);
2341     C2F(drot)(&iFour, &pdblTemp[0][0], &iThree, &pdblTemp[1][0], &iThree, &dblD, &dblE);
2342
2343     //z1
2344     dgivs(pdblTemp[1][3], pdblTemp[1][2], &dblD, &dblE);
2345     C2F(drot)(&iTwo, &pdblTemp[0][3], &iOne, &pdblTemp[0][2], &iOne, &dblD, &dblE);
2346     C2F(drot)(&iIndex3, &_pdblIn[(iIndex3 - 1) * _iLeadDim], &iOne, &_pdblIn[(iIndex2 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2347     C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex3 - 1) * _iLeadDim], &iOne, &_pdblOut[(iIndex2 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2348
2349     //z1'
2350     C2F(drot)(&dblTemp2, &_pdblIn[(iIndex3 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim], &iOne, &dblD, &dblE);
2351
2352     //z2
2353     dgivs(pdblTemp[1][3], pdblTemp[1][1], &dblD, &dblE);
2354     C2F(drot)(&iTwo, &pdblTemp[0][3], &iOne, &pdblTemp[0][1], &iOne, &dblD, &dblE);
2355     C2F(drot)(&iIndex3, &_pdblIn[(iIndex3 - 1) * _iLeadDim], &iOne, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2356     C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex3 - 1) * _iLeadDim], &iOne, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2357
2358     //z2'
2359     pdblTemp[1][3] *=  dblD;
2360     C2F(drot)(&dblTemp2, &_pdblIn[(iIndex3 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &_iMax, &dblD, &dblE);
2361
2362     //z3
2363     dgivs(pdblTemp[0][2], pdblTemp[0][1], &dblD, &dblE);
2364     C2F(drot)(&iOne, &pdblTemp[0][2], &iOne, &pdblTemp[0][1], &iOne, &dblD, &dblE);
2365     C2F(drot)(&iIndex3, &_pdblIn[(iIndex2 - 1) * _iLeadDim], &iOne, &_pdblIn[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2366     C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex2 - 1) * _iLeadDim], &iOne, &_pdblOut[(iIndex1 - 1) * _iLeadDim], &iOne, &dblD, &dblE);
2367
2368     //z3'
2369     pdblTemp[1][3] *=  dblD;
2370     C2F(drot)(&dblTemp2, &_pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[(iIndex1 - 1) + _iPos * _iLeadDim], &_iMax, &dblD, &dblE);
2371
2372     //z4
2373     dgivs(pdblTemp[0][2], pdblTemp[0][0], &dblD, &dblE);
2374     C2F(drot)(&iIndex3, &_pdblIn[(iIndex2 - 1) * _iLeadDim], &iOne, &_pdblIn[_iPos * _iLeadDim], &iOne, &dblD, &dblE);
2375     C2F(drot)(&_iLeadDim, &_pdblOut[(iIndex2 - 1) * _iLeadDim], &iOne, &_pdblOut[_iPos * _iLeadDim], &iOne, &dblD, &dblE);
2376
2377     //z4'
2378     C2F(drot)(&dblTemp2, &_pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim], &_iMax, &_pdblIn[_iPos + _iPos * _iLeadDim], &_iMax, &dblD, &dblE);
2379
2380     //zeroes negligible elements
2381     _pdblIn[(iIndex2 - 1) + _iPos * _iLeadDim]                  = 0;
2382     _pdblIn[(iIndex3 - 1) + _iPos * _iLeadDim]                  = 0;
2383     _pdblIn[(iIndex2 - 1) + (iIndex1 - 1) * _iLeadDim]  = 0;
2384     _pdblIn[(iIndex3 - 1) + (iIndex1 - 1) * _iLeadDim]  = 0;
2385
2386     return 0;
2387 }
2388
2389 /*
2390         purpose
2391                  this routine constructs the givens transformation
2392
2393                                    ( sc  ss )
2394                            g = (        ),   sc**2+ss**2 = 1. ,
2395                                    (-ss  sc )
2396
2397                  which zeros the second entry of the 2-vector (sa,sb)**t
2398                  this routine is a modification of the blas routine srotg
2399                  (algorithm 539) in order to leave the arguments sa and sb
2400                  unchanged
2401
2402         syntax
2403
2404                 subroutine giv(sa,sb,sc,ss)
2405                 double precision sa,sb,sc,ss
2406         auxiliary routines
2407                 sqrt abs (fortran)
2408 */
2409 int dgivs(double _dblA, double _dblB, double *_pdblSC, double *_pdblSS)
2410 {
2411     double dblU = 0;
2412     double dblR = 0;
2413     double dblV = 0;
2414     if (dabss(_dblA) > dabss(_dblB))
2415     {
2416         dblU            = _dblA + _dblA;
2417         dblV            = _dblB / dblU;
2418         dblR            = dsqrts(0.25 + dblV * dblV) * dblU;
2419         *_pdblSC        = _dblA / dblR;
2420         *_pdblSS        = dblV * (*_pdblSC + *_pdblSC);
2421     }
2422     else if (_dblB != 0)
2423     {
2424         dblU            = _dblB + _dblB;
2425         dblV            = _dblA / dblU;
2426         dblR            = dsqrts(0.25 + dblV * dblV) * dblU;
2427         *_pdblSS        = _dblB / dblR;
2428         *_pdblSC        = dblV * (*_pdblSS + *_pdblSS);
2429     }
2430     else
2431     {
2432         //dblB == 0
2433         *_pdblSC        = 0;
2434         *_pdblSS        = 0;
2435     }
2436     return 0;
2437 }
2438
2439 /*
2440 purpose
2441
2442     given the upper hessenberg matrix a with a 2x2 block
2443     starting at a(l,l), split determines if the
2444     corresponding eigenvalues are real or complex, if they
2445     are real, a rotation is determined that reduces the
2446     block to upper triangular form with the eigenvalue
2447     of largest absolute value appearing first.  the
2448     rotation is accumulated in v. the eigenvalues (real
2449     or complex) are returned in e1 and e2.
2450 syntax
2451
2452     subroutine split(a, v, n, l, e1, e2, na, nv)
2453
2454     double precision a,v,e1,e2
2455     integer n,l,na,nv
2456     dimension a(na,n),v(nv,n)
2457
2458     starred parameters are  altered by the subroutine
2459
2460        *a        the upper hessenberg matrix whose 2x2
2461                  block is to be dsplit.
2462        *v        the array in which the dsplitting trans-
2463                  formation is to be accumulated.
2464         n        the order of the matrix a.
2465         l        the position of the 2x2 block.
2466        *e1       on return if the eigenvalues are complex
2467        *e2       e1 contains their common real part and
2468                  e2 contains the positive imaginary part.
2469                  if the eigenvalues are real. e1 contains
2470                  the one largest in absolute value and f2
2471                  contains the other one.
2472        na        the first dimension of the array a.
2473        nv        the first dimension of the array v.
2474 auxiliary routines
2475     abs sqrt (fortran)
2476 */
2477
2478 int dsplits(double *_pdblVal, double *_pdblSplit, int _iSize, int _iPos, double *_pdblE1, double *_pdblE2, int _iLeadDimVal, int _iLeadDimSplit)
2479 {
2480     double dblZero      = 0;
2481     double dblTwo       = 2;
2482     double dblP = 0, dblQ = 0, dblR = 0, dblS = 0, dblT = 0, dblU = 0;
2483     double dblW = 0, dblX = 0, dblY = 0, dblZ = 0;
2484
2485     int iLoop1  = 0;
2486     int iLoop2  = 0;
2487     int iIndex  = 0;
2488
2489     iIndex      = _iPos + 1;
2490
2491     dblX        =       _pdblVal[(iIndex - 1) + (iIndex - 1) * _iLeadDimVal];
2492     dblY        =       _pdblVal[_iPos + _iPos * _iLeadDimVal];
2493     dblW        =       _pdblVal[_iPos + (iIndex - 1) * _iLeadDimVal] *
2494                 _pdblVal[(iIndex - 1) + _iPos * _iLeadDimVal];
2495     dblP        = (dblY / dblX) / dblTwo;
2496     dblQ        = dblP * dblP + dblW;
2497
2498     if (dblQ < 0)
2499     {
2500         //complex eigenvalue.
2501         *_pdblE1        = dblP + dblX;
2502         *_pdblE2        = dsqrts(-dblQ);
2503         return 0;
2504     }
2505
2506     //L10
2507     //two real eigenvalues. set up transformation
2508     dblZ        = dsqrts(dblQ);
2509     if (dblP < dblZero)
2510     {
2511         dblZ    = dblP + dblZ;
2512     }
2513     else
2514     {
2515         dblZ    = dblP - dblZ;
2516     }
2517
2518     if (dblZ == dblZero)
2519     {
2520         dblR    = dblZero;
2521     }
2522     else
2523     {
2524         dblR    = -dblW / dblZ;
2525     }
2526
2527     if (dabss(dblX + dblZ) >= dabss(dblX + dblR))
2528     {
2529         dblZ    = dblR;
2530     }
2531
2532     dblY        = dblY - dblX - dblZ;
2533     dblX        = -dblZ;
2534     dblT        = _pdblVal[_iPos + (iIndex - 1) * _iLeadDimVal];
2535     dblU        = _pdblVal[(iIndex - 1) + _iPos * _iLeadDimVal];
2536
2537     if (dabss(dblY) + dabss(dblU) <= dabss(dblT) + dabss(dblX))
2538     {
2539         dblQ    = dblX;
2540         dblP    = dblT;
2541     }
2542     else
2543     {
2544         dblQ    = dblU;
2545         dblP    = dblY;
2546     }
2547
2548     dblR        = dsqrts(dblP * dblP + dblQ * dblQ);
2549     if (dblR <= dblZero)
2550     {
2551         *_pdblE1        = _pdblVal[_iPos + _iPos * _iLeadDimVal];
2552         *_pdblE2        = _pdblVal[(iIndex - 1) + (iIndex - 1) * _iLeadDimVal];
2553         _pdblVal[(iIndex - 1) + _iPos * _iLeadDimVal]   = dblZero;
2554         return 0;
2555     }
2556
2557     dblP        /= dblR;
2558     dblQ        /= dblR;
2559
2560     //premultiply
2561     for (iLoop1 = _iPos ; iLoop1 < _iSize ; iLoop1++)
2562     {
2563         dblZ    = _pdblVal[_iPos + iLoop1 * _iLeadDimVal];
2564         _pdblVal[_iPos + iLoop1 * _iLeadDimVal]                 = dblP * dblZ +
2565                 dblQ * _pdblVal[(iIndex - 1) + iLoop1 * _iLeadDimVal];
2566         _pdblVal[(iIndex - 1) + iLoop1 * _iLeadDimVal]  = dblP *
2567                 _pdblVal[(iIndex - 1) + iLoop1 * _iLeadDimVal] - dblQ * dblZ;
2568     }
2569
2570     //postmultiply
2571     for (iLoop1 = 0 ; iLoop1 < iIndex ; iLoop1++)
2572     {
2573         dblZ    = _pdblVal[iLoop1 + _iPos * _iLeadDimVal];
2574         _pdblVal[iLoop1 + _iPos * _iLeadDimVal]                 = dblP * dblZ +
2575                 dblQ * _pdblVal[iLoop1 + (iIndex - 1) * _iLeadDimVal];
2576         _pdblVal[iLoop1 + (iIndex - 1) * _iLeadDimVal]  = dblP *
2577                 _pdblVal[iLoop1 + (iIndex - 1) * _iLeadDimVal] - dblQ * dblZ;
2578     }
2579
2580     //accumulate the transformation in v
2581     for (iLoop1 = 0 ; iLoop1 < _iSize ; iLoop1++)
2582     {
2583         dblZ    = _pdblSplit[iLoop1 + _iPos * _iLeadDimVal];
2584         _pdblSplit[iLoop1 + _iPos * _iLeadDimVal]                       = dblP * dblZ +
2585                 dblQ * _pdblSplit[iLoop1 + (iIndex - 1) * _iLeadDimVal];
2586         _pdblSplit[iLoop1 + (iIndex - 1) * _iLeadDimVal]        = dblP *
2587                 _pdblSplit[iLoop1 + (iIndex - 1) * _iLeadDimVal] - dblQ * dblZ;
2588     }
2589
2590     _pdblVal[(iIndex - 1) + _iPos * _iLeadDimVal]       = dblZero;
2591     *_pdblE1    = _pdblVal[_iPos + _iPos * _iLeadDimVal];
2592     *_pdblE2    = _pdblVal[(iIndex - 1) + (iIndex - 1) * _iLeadDimVal];
2593     return 0;
2594 }
2595
2596
2597 /*
2598 C2F(dgemm)("n", "n", &iRows1, &iCols2, &iCols1, &dblOne, pReal1 , &iRows1 , pReal2, &iRows2, &dblZero, pReturnReal ,&iRows1);
2599 if(iComplex1 == 1)
2600         C2F(dgemm)("n", "n", &iRows1, &iCols2, &iCols1, &dblOne, pImg1 , &iRows1 , pReal2, &iRows2, &dblZero, pReturnImg ,&iRows1);
2601 if(iComplex2 == 1)
2602         C2F(dgemm)("n", "n", &iRows1, &iCols2, &iCols1, &dblOne, pReal1 , &iRows1 , pImg2, &iRows2, &dblZero, pReturnImg ,&iRows1);
2603 */
2604
2605 int dexpms2(double *_pdblReal, double *_pdblReturnReal, int _iLeadDim)
2606 {
2607     int iRet = zexpms2(_pdblReal, NULL, _pdblReturnReal, NULL, _iLeadDim);
2608     return iRet;
2609 }
2610
2611 int zexpms2(double *_pdblReal, double *_pdblImg, double *_pdblReturnReal, double *_pdblReturnImg, int _iLeadDim)
2612 {
2613     double dblRcond = 0;
2614     int iRet                            = 0;
2615     int iIndex1                 = 0;
2616     int iMax                            = 0;
2617     int iFlag                           = 0;
2618     int iLoop1                  = 0;
2619     int iSquare                 = 0;
2620     int iOne                            = 1;
2621
2622     int iComplex = 0;
2623
2624     double dblZero      = 0;
2625     double dblOne       = 1;
2626
2627     double dblExp       = 0;
2628     double dblS         = 0;
2629     double dblS2        = 0;
2630     double dblCst       = 0.5;
2631
2632     double *pdblMatrixRealA             = NULL;//A'
2633     double *pdblMatrixRealX             = NULL;//X
2634     double *pdblMatrixRealD             = NULL;//D
2635     double *pdblMatrixRealcX    = NULL;//cX
2636     double *pdblMatrixRealcA    = NULL;//cX
2637     double *pdblMatrixRealEye   = NULL;//Eye
2638     double *pdblMatrixRealTemp  = NULL;//Temp
2639     double *pdblMatrixRealTemp2 = NULL;//Temp2
2640
2641     double *pdblMatrixImgA              = NULL;//A'
2642     double *pdblMatrixImgX              = NULL;//X
2643     double *pdblMatrixImgD              = NULL;//D
2644     double *pdblMatrixImgcX             = NULL;//cX
2645     double *pdblMatrixImgcA             = NULL;//cX
2646     double *pdblMatrixImgEye    = NULL;//Eye
2647     double *pdblMatrixImgTemp   = NULL;//Temp
2648     double *pdblMatrixImgTemp2  = NULL;//Temp2
2649
2650     if (_pdblImg != NULL)
2651     {
2652         iComplex = 1;
2653     }
2654
2655     iSquare = _iLeadDim * _iLeadDim;
2656
2657     pdblMatrixRealA                     = (double*)malloc(sizeof(double) * iSquare);
2658     pdblMatrixRealX                     = (double*)malloc(sizeof(double) * iSquare);
2659     pdblMatrixRealD                     = (double*)malloc(sizeof(double) * iSquare);
2660     pdblMatrixRealcX            = (double*)malloc(sizeof(double) * iSquare);
2661     pdblMatrixRealcA            = (double*)malloc(sizeof(double) * iSquare);
2662     pdblMatrixRealEye           = (double*)malloc(sizeof(double) * iSquare);
2663     pdblMatrixRealTemp          = (double*)malloc(sizeof(double) * iSquare);
2664     pdblMatrixRealTemp2         = (double*)malloc(sizeof(double) * iSquare);
2665
2666     if (iComplex)
2667     {
2668         pdblMatrixImgA                  = (double*)malloc(sizeof(double) * iSquare);
2669         pdblMatrixImgX                  = (double*)malloc(sizeof(double) * iSquare);
2670         pdblMatrixImgD                  = (double*)malloc(sizeof(double) * iSquare);
2671         pdblMatrixImgcX                 = (double*)malloc(sizeof(double) * iSquare);
2672         pdblMatrixImgcA                 = (double*)malloc(sizeof(double) * iSquare);
2673         pdblMatrixImgEye                = (double*)malloc(sizeof(double) * iSquare);
2674         pdblMatrixImgTemp               = (double*)malloc(sizeof(double) * iSquare);
2675         pdblMatrixImgTemp2              = (double*)malloc(sizeof(double) * iSquare);
2676
2677         memset(pdblMatrixImgEye, 0x00, sizeof(double) * iSquare);
2678     }
2679
2680
2681     // Scale A by power of 2 so that its norm is < 1/2 .
2682     dfrexps(dblGetMatrixInfiniteNorm(_pdblReal, _pdblImg, _iLeadDim, _iLeadDim), &dblExp);
2683     dblS        = Max(0, dblExp + 1);
2684     dblS2       = pow(2, dblS);
2685
2686     //A = A./2^s
2687     if (iComplex)
2688     {
2689         iRightDivisionComplexMatrixByRealMatrix(
2690             _pdblReal,                  _pdblImg        ,               1,
2691             &dblS2,                             0,
2692             pdblMatrixRealA,    pdblMatrixImgA, 1,      iSquare);
2693     }
2694     else
2695     {
2696         iRightDivisionRealMatrixByRealMatrix(
2697             _pdblReal,                  1,
2698             &dblS2,                             0,
2699             pdblMatrixRealA,    1,      iSquare);
2700     }
2701
2702     // Pade approximation for exp(A)
2703     //X = A
2704     C2F(dcopy)(&iSquare, pdblMatrixRealA, &iOne, pdblMatrixRealX, &iOne );
2705     if (iComplex)
2706     {
2707         C2F(dcopy)(&iSquare, pdblMatrixImgA, &iOne, pdblMatrixImgX, &iOne );
2708     }
2709
2710
2711     deyes(pdblMatrixRealEye, _iLeadDim, _iLeadDim);
2712
2713
2714     //cA = A * c
2715     if (iComplex)
2716         iMultiRealScalarByComplexMatrix(
2717             dblCst,
2718             pdblMatrixRealA, pdblMatrixImgA, _iLeadDim, _iLeadDim,
2719             pdblMatrixRealcA, pdblMatrixImgcA);
2720     else
2721         iMultiRealScalarByRealMatrix(
2722             dblCst,
2723             pdblMatrixRealA, _iLeadDim, _iLeadDim,
2724             pdblMatrixRealcA);
2725     /*
2726         C2F(dcopy)(&iSquare, pdblMatrixRealA, &iOne, pdblMatrixRealcA, &iOne);
2727         C2F(dscal)(&iSquare, &dblCst, pdblMatrixRealcA, &iOne);
2728     */
2729
2730     //E = Eye + cA
2731     vDadd(iSquare, pdblMatrixRealEye, pdblMatrixRealcA, 1, 1, _pdblReturnReal);
2732     if (iComplex)
2733     {
2734         vDadd(iSquare, pdblMatrixImgEye, pdblMatrixImgcA, 1, 1, _pdblReturnImg);
2735     }
2736
2737     //D = Eye - cA
2738     vDless(iSquare, pdblMatrixRealEye, pdblMatrixRealcA, 1, 1, pdblMatrixRealD);
2739     if (iComplex)
2740     {
2741         vDless(iSquare, pdblMatrixImgEye, pdblMatrixImgcA, 1, 1, pdblMatrixImgD);
2742     }
2743
2744     iMax        = 6;
2745     iFlag       = 1;
2746
2747     for (iLoop1 = 2 ; iLoop1 <= iMax ; iLoop1++)
2748     {
2749         dblCst  = dblCst * (iMax - iLoop1 + 1 ) / (iLoop1 * (2 * iMax - iLoop1 + 1));
2750
2751         //Temp = X
2752         C2F(dcopy)(&iSquare, pdblMatrixRealX, &iOne, pdblMatrixRealTemp, &iOne);
2753         if (iComplex)
2754         {
2755             C2F(dcopy)(&iSquare, pdblMatrixImgX, &iOne, pdblMatrixImgTemp, &iOne);
2756         }
2757         /*              for(iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++)
2758                                 pdblMatrixRealTemp[iIndex1]     = pdblMatrixRealX[iIndex1];
2759         */
2760         //X = A * Temp;
2761         if (iComplex)
2762         {
2763             iMultiComplexMatrixByComplexMatrix(
2764                 pdblMatrixRealA, pdblMatrixImgA, _iLeadDim, _iLeadDim,
2765                 pdblMatrixRealTemp, pdblMatrixImgTemp, _iLeadDim, _iLeadDim,
2766                 pdblMatrixRealX, pdblMatrixImgX);
2767         }
2768         else
2769         {
2770             iMultiRealMatrixByRealMatrix(
2771                 pdblMatrixRealA, _iLeadDim, _iLeadDim,
2772                 pdblMatrixRealTemp, _iLeadDim, _iLeadDim,
2773                 pdblMatrixRealX);
2774         }
2775         /*              C2F(dgemm)("n", "n", &_iLeadDim, &_iLeadDim, &_iLeadDim, &dblOne,
2776                                 pdblMatrixRealA, &_iLeadDim ,
2777                                 pdblMatrixRealTemp, &_iLeadDim, &dblZero,
2778                                 pdblMatrixRealX ,&_iLeadDim);
2779         */
2780         //cX = c * X
2781         if (iComplex)
2782             iMultiRealScalarByComplexMatrix(
2783                 dblCst,
2784                 pdblMatrixRealX, pdblMatrixImgX, _iLeadDim, _iLeadDim,
2785                 pdblMatrixRealcX, pdblMatrixImgcX);
2786         else
2787             iMultiRealScalarByRealMatrix(
2788                 dblCst,
2789                 pdblMatrixRealX, _iLeadDim, _iLeadDim,
2790                 pdblMatrixRealcX);
2791         /*              C2F(dcopy)(&iSquare, pdblMatrixRealX, &iOne, pdblMatrixRealcX, &iOne);
2792                         C2F(dscal)(&iSquare, &dblCst, pdblMatrixRealcX, &iOne);
2793         */
2794
2795         //E = E + cX
2796         vDadd(iSquare, _pdblReturnReal, pdblMatrixRealcX, 1, 1, _pdblReturnReal);
2797         if (iComplex)
2798         {
2799             vDadd(iSquare, _pdblReturnImg, pdblMatrixImgcX, 1, 1, _pdblReturnImg);
2800         }
2801
2802
2803         if (iFlag == 1) //D = D + cX
2804         {
2805             vDadd(iSquare, pdblMatrixRealD, pdblMatrixRealcX, 1, 1, pdblMatrixRealD);
2806             if (iComplex)
2807             {
2808                 vDadd(iSquare, pdblMatrixImgD, pdblMatrixImgcX, 1, 1, pdblMatrixImgD);
2809             }
2810         }
2811         else //D = D - cX
2812         {
2813             vDless(iSquare, pdblMatrixRealD, pdblMatrixRealcX, 1, 1, pdblMatrixRealD);
2814             if (iComplex)
2815             {
2816                 vDless(iSquare, pdblMatrixImgD, pdblMatrixImgcX, 1, 1, pdblMatrixImgD);
2817             }
2818         }
2819
2820         //Toggle iFlag
2821         iFlag = !iFlag;
2822     }
2823
2824     //Temp = E
2825     C2F(dcopy)(&iSquare, _pdblReturnReal, &iOne, pdblMatrixRealTemp, &iOne);
2826     if (iComplex)
2827     {
2828         C2F(dcopy)(&iSquare, _pdblReturnImg, &iOne, pdblMatrixImgTemp, &iOne);
2829         /*      for(iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++)
2830                         pdblMatrixRealTemp[iIndex1]     = _pdblReturnReal[iIndex1];
2831         */
2832     }
2833
2834     //E = D\E
2835     if (iComplex)
2836     {
2837         iRet = iLeftDivisionOfComplexMatrix(
2838                    pdblMatrixRealD,     pdblMatrixImgD,         _iLeadDim, _iLeadDim,
2839                    pdblMatrixRealTemp,  pdblMatrixImgTemp,      _iLeadDim, _iLeadDim,
2840                    _pdblReturnReal,     _pdblReturnImg,         _iLeadDim, _iLeadDim, &dblRcond);
2841     }
2842     else
2843     {
2844         iRet = iLeftDivisionOfRealMatrix(
2845                    pdblMatrixRealD,     _iLeadDim, _iLeadDim,
2846                    pdblMatrixRealTemp,  _iLeadDim, _iLeadDim,
2847                    _pdblReturnReal,     _iLeadDim, _iLeadDim, &dblRcond);
2848     }
2849
2850     if (iRet < 0)
2851     {
2852         switch (iRet)
2853         {
2854             case -1 :
2855                 //Scilab 6
2856                 //sprintf(C2F(cha1).buf, "%1.4E", dblRcond);
2857                 //Msgs(5,1);
2858                 break;
2859             case -2 :
2860                 //Scilab 6
2861                 //Msgs(9, (int)dblRcond);
2862                 break;
2863             default :
2864                 break;
2865         }
2866     }
2867     // Undo scaling by repeated squaring
2868     for (iLoop1 = 0 ; iLoop1 < dblS ; iLoop1++)
2869     {
2870         //Temp = E
2871         C2F(dcopy)(&iSquare, _pdblReturnReal, &iOne, pdblMatrixRealTemp, &iOne);
2872         //Temp2 = E
2873         C2F(dcopy)(&iSquare, _pdblReturnReal, &iOne, pdblMatrixRealTemp2, &iOne);
2874         if (iComplex)
2875         {
2876             //Temp = E
2877             C2F(dcopy)(&iSquare, _pdblReturnImg, &iOne, pdblMatrixImgTemp, &iOne);
2878             //Temp2 = E
2879             C2F(dcopy)(&iSquare, _pdblReturnImg, &iOne, pdblMatrixImgTemp2, &iOne);
2880         }
2881
2882
2883         /*              for(iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++)
2884                                 pdblMatrixRealTemp[iIndex1]             = _pdblReturnReal[iIndex1];
2885         */
2886         /*              for(iIndex1 = 0 ; iIndex1 < iSquare ; iIndex1++)
2887                                 pdblMatrixRealTemp2[iIndex1]    = _pdblReturnReal[iIndex1];
2888         */
2889         // E = E*E
2890         if (iComplex)
2891             iMultiComplexMatrixByComplexMatrix(
2892                 pdblMatrixRealTemp,             pdblMatrixImgTemp,      _iLeadDim, _iLeadDim,
2893                 pdblMatrixRealTemp2,    pdblMatrixImgTemp2,     _iLeadDim, _iLeadDim,
2894                 _pdblReturnReal,                _pdblReturnImg);
2895         else
2896             iMultiRealMatrixByRealMatrix(
2897                 pdblMatrixRealTemp,             _iLeadDim, _iLeadDim,
2898                 pdblMatrixRealTemp2,    _iLeadDim, _iLeadDim,
2899                 _pdblReturnReal);
2900         /*
2901                         C2F(dgemm)("n", "n", &_iLeadDim, &_iLeadDim, &_iLeadDim, &dblOne,
2902                                 pdblMatrixRealTemp,             &_iLeadDim ,
2903                                 pdblMatrixRealTemp2,    &_iLeadDim, &dblZero,
2904                                 _pdblReturnReal ,               &_iLeadDim);
2905         */
2906     }
2907
2908     free(pdblMatrixRealA);
2909     free(pdblMatrixRealX);
2910     free(pdblMatrixRealD);
2911     free(pdblMatrixRealcX);
2912     free(pdblMatrixRealcA);
2913     free(pdblMatrixRealEye);
2914     free(pdblMatrixRealTemp);
2915     free(pdblMatrixRealTemp2);
2916
2917     if (iComplex)
2918     {
2919         free(pdblMatrixImgA);
2920         free(pdblMatrixImgX);
2921         free(pdblMatrixImgD);
2922         free(pdblMatrixImgcX);
2923         free(pdblMatrixImgcA);
2924         free(pdblMatrixImgEye);
2925         free(pdblMatrixImgTemp);
2926         free(pdblMatrixImgTemp2);
2927     }
2928
2929     return 0;
2930 }
2931
2932 double dblGetMatrixInfiniteNorm(double *_pdblReal, double *_pdblImg, int _iRows, int _iCols)
2933 {
2934     int iIndex1 = 0, iIndex2 = 0;
2935     double dblTemp = 0;
2936     double dblRef = 0;
2937
2938     if (_pdblImg == NULL)
2939     {
2940         for (iIndex1 = 0 ; iIndex1 < _iRows ; iIndex1++)
2941         {
2942             dblTemp = 0;
2943             for (iIndex2 = 0 ; iIndex2 < _iCols ; iIndex2++)
2944             {
2945                 dblTemp += _pdblReal[iIndex1 + iIndex2 * _iRows];
2946             }
2947             if (dblTemp > dblRef)
2948             {
2949                 dblRef = dblTemp;
2950             }
2951         }
2952     }
2953     else
2954     {
2955         for (iIndex1 = 0 ; iIndex1 < _iRows ; iIndex1++)
2956         {
2957             dblTemp = 0;
2958             for (iIndex2 = 0 ; iIndex2 < _iCols ; iIndex2++)
2959             {
2960                 dblTemp += dpythags(_pdblReal[iIndex1 + iIndex2 * _iRows], _pdblImg[iIndex1 + iIndex2 * _iRows]);
2961             }
2962             if (dblTemp > dblRef)
2963             {
2964                 dblRef = dblTemp;
2965             }
2966         }
2967     }
2968     return dblRef;
2969 }