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