Merge remote-tracking branch 'origin/master' into jit
[scilab.git] / scilab / modules / ast / src / c / operations / matrix_power.c
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2009 - DIGITEO - Antoine ELIAS
4 *
5 *  This file must be used under the terms of the CeCILL.
6 *  This source file is licensed as described in the file COPYING, which
7 *  you should have received as part of this distribution.  The terms
8 *  are also available at
9 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <string.h>
16 #include "matrix_power.h"
17 #include "matrix_multiplication.h"
18 //#include "unsfdcopy.h"
19
20 #include "elementary_functions.h"
21 #include "elem_common.h"
22 #include "invert_matrix.h"
23 #include "sciprint.h"
24 #include "localization.h"
25 #include "configvariable_interface.h"
26 #include "Sciwarning.h"
27
28 /*
29 r : real part
30 c : imaginary part
31 E : element of
32 R : real
33 Z : integer
34 C : complex
35 * : less 0
36 + : >= 0
37 - : <= 0
38 */
39 extern int C2F(zgemm) (char *_pstTransA, char *_pstTransB, int *_piN, int *_piM, int *_piK, double *_pdblAlpha, double *_pdblA, int *_piLdA,
40                        double *_pdblB, int *_piLdB, double *_pdblBeta, double *_pdblC, int *_piLdC);
41
42 /*ddpowe*/
43 int iPowerRealScalarByRealScalar(
44     double _dblReal1,
45     double _dblReal2,
46     double *_pdblRealOut, double *_pdblImgOut, int *_piComplex)
47 {
48     //exposant is an integer
49     if ((int)_dblReal2 == _dblReal2)
50     {
51         //dipowe
52         int iReal2 = (int)_dblReal2;
53
54         if (iReal2 == 1)
55         {
56             //R ^ 1
57             *_pdblRealOut = _dblReal1;
58             *_pdblImgOut        = 0;
59             *_piComplex         = 0;
60         }
61         else if (iReal2 == 0)
62         {
63             //R ^ 0
64             *_pdblRealOut = 1;
65             *_pdblImgOut        = 0;
66             *_piComplex         = 0;
67         }
68         else if (iReal2 < 0)
69         {
70             //R ^ Z-
71             if (_dblReal1 != 0)
72             {
73                 //R* ^ Z-
74                 *_pdblRealOut = pow(_dblReal1, iReal2);
75                 *_pdblImgOut    = 0;
76                 *_piComplex             = 0;
77             }
78             else
79             {
80                 //0 ^ 0
81                 //FIXME : ieee
82                 //generate +Inf
83                 double dblZero  = 0.0;
84                 *_pdblRealOut           = 1.0 / (dblZero);
85                 *_pdblImgOut            = 0;
86                 *_piComplex                     = 0;
87             }
88         }
89         else
90         {
91             //R ^ Z*+ ( N* )
92             *_pdblRealOut = pow(_dblReal1, iReal2);
93             *_pdblImgOut        = 0;
94             *_piComplex         = 0;
95         }
96     }
97     else
98     {
99         if (_dblReal1 > 0)
100         {
101             //R*+ ^ R
102             *_pdblRealOut = pow(_dblReal1, _dblReal2);
103             *_pdblImgOut        = 0;
104             *_piComplex         = 0;
105         }
106         else if (_dblReal1 < 0)
107         {
108             //R*- ^ R
109             double dblRealTemp  = 0;
110             double dblImgTemp           = 0;
111
112             wlog(_dblReal1, 0, &dblRealTemp, &dblImgTemp);
113             dblRealTemp                                 = dexps(dblRealTemp * _dblReal2);
114             dblImgTemp                                  = dblImgTemp * _dblReal2;
115             *_pdblRealOut                               = dblRealTemp * dcoss(dblImgTemp);
116             *_pdblImgOut                                = dblRealTemp * dsins(dblImgTemp);
117             *_piComplex                                 = 1;
118         }
119         else if (_dblReal1 == 0)
120         {
121             //0 ^ R
122             if (_dblReal2 < 0)
123             {
124                 //0 ^ R*-
125                 //FIXME : ieee
126                 //generate +Inf
127                 double dblZero  = 0.0;
128                 *_pdblRealOut           = 1.0 / (dblZero);
129                 *_pdblImgOut            = 0;
130                 *_piComplex                     = 0;
131             }
132             else if (_dblReal2 == 0)
133             {
134                 //0 ^ 0
135                 //never call, pass in R ** 0 before
136                 *_pdblRealOut           = 1;
137                 *_pdblImgOut            = 0;
138                 *_piComplex                     = 0;
139             }
140             else if (_dblReal2 > 0)
141             {
142                 //0 ^ R*+
143                 *_pdblRealOut           = 0;
144                 *_pdblImgOut            = 0;
145                 *_piComplex                     = 0;
146             }
147             else //_dblReal2 is NaN
148             {
149                 *_pdblRealOut           = _dblReal2;
150                 *_pdblImgOut            = 0;
151                 *_piComplex                     = 0;
152             }
153         }
154         else //_dblReal1 is NaN
155         {
156             *_pdblRealOut               = _dblReal1;
157             *_pdblImgOut                = 0;
158             *_piComplex                 = 0;
159         }
160     }
161     return 0;
162 }
163
164 /*dwpowe*/
165 int iPowerRealScalarByComplexScalar(
166     double _dblReal1,
167     double _dblReal2, double _dblImg2,
168     double* _pdblRealOut, double* _pdblImgOut)
169 {
170     if (_dblImg2 == 0)
171     {
172         //R ^ R
173         int iComplex = 0;
174         iPowerRealScalarByRealScalar(
175             _dblReal1,
176             _dblReal2,
177             _pdblRealOut, _pdblImgOut, &iComplex);
178     }
179     else
180     {
181         //R ^ C
182         if (_dblReal1 != 0)
183         {
184             //R* ^ C
185             double dblRealTemp  = 0;
186             double dblImgTemp           = 0;
187
188             wlog(_dblReal1, 0, &dblRealTemp, &dblImgTemp);
189             C2F(wmul)(&dblRealTemp, &dblImgTemp, &_dblReal2, &_dblImg2, &dblRealTemp, &dblImgTemp);
190             dblRealTemp                                 = dexps(dblRealTemp);
191             *_pdblRealOut                               = dblRealTemp * dcoss(dblImgTemp);
192             *_pdblImgOut                                = dblRealTemp * dsins(dblImgTemp);
193         }
194         else
195         {
196             //0 ^ C
197             if (_dblReal2 > 0)
198             {
199                 //0 ^ (r E R*+ ) & ( c E R )
200                 *_pdblRealOut                           = 0;
201                 *_pdblImgOut                            = 0;
202             }
203             else if (_dblReal2 < 0)
204             {
205                 //0 ^ (r E R*- ) & ( c E R )
206                 //FIXME : ieee
207                 //generate +Inf
208                 double dblZero  = 0.0;
209                 *_pdblRealOut           = 1.0 / (dblZero);
210                 *_pdblImgOut            = 0;
211             }
212             else //_dblReal2 == 0, NaN
213             {
214                 //0 ^ (r = 0 ) & ( c E R )
215                 *_pdblRealOut                           = 1;
216                 *_pdblImgOut                            = 0;
217             }
218         }
219     }
220     return 0;
221 }
222
223 /*wdpowe*/
224 int iPowerComplexScalarByRealScalar(
225     double _dblReal1, double _dblImg1,
226     double _dblReal2,
227     double* _pdblRealOut, double* _pdblImgOut)
228 {
229     if ((int)_dblReal2 == _dblReal2)
230     {
231         //C ^ Z
232         if (_dblReal2 == 0)
233         {
234             //C ^ 0
235             *_pdblRealOut                               = 1;
236             *_pdblImgOut                                = 0;
237         }
238         else if (_dblReal2 < 0)
239         {
240             //C ^ Z*-
241             if (dabss(_dblReal1) + dabss(_dblImg1) != 0) //_dblReal1 != 0 || _dblImg1 != 0 ?
242             {
243                 int i = 0;
244                 double dblOne                           = 1;
245                 double dblZero                  = 0;
246                 double dblRealTemp      = 0;
247                 double dblImgTemp               = 0;
248
249                 C2F(wdiv)(&dblOne, &dblZero, &_dblReal1, &_dblImg1, _pdblRealOut, _pdblImgOut);
250
251                 dblRealTemp                                     = *_pdblRealOut;
252                 dblImgTemp                                      = *_pdblImgOut;
253
254                 for (i = 2 ; i <= dabss(_dblReal2) ; i++)
255                 {
256                     C2F(wmul)(&dblRealTemp, &dblImgTemp, _pdblRealOut, _pdblImgOut, _pdblRealOut, _pdblImgOut);
257                 }
258             }
259             else
260             {
261                 //FIXME : ieee
262                 //generate +Inf
263                 double dblZero  = 0.0;
264                 *_pdblRealOut           = 1.0 / (dblZero);
265                 *_pdblImgOut            = 0;
266             }
267         }
268         else
269         {
270             //C ^ Z*+
271             int i                                                               = 0;
272             double dblRealTemp  = 0;
273             double dblImgTemp           = 0;
274
275             *_pdblRealOut                       = _dblReal1;
276             *_pdblImgOut                        = _dblImg1;
277             dblRealTemp                         = *_pdblRealOut;
278             dblImgTemp                          = *_pdblImgOut;
279
280             for (i = 2 ; i <= dabss(_dblReal2) ; i++)
281             {
282                 C2F(wmul)(&dblRealTemp, &dblImgTemp, _pdblRealOut, _pdblImgOut, _pdblRealOut, _pdblImgOut);
283             }
284         }
285     }
286     else
287     {
288         if (dabss(_dblReal1) + dabss(_dblImg1) != 0)
289         {
290             //C ^ R
291             double dblRealTemp  = 0;
292             double dblImgTemp           = 0;
293
294             wlog(_dblReal1, _dblImg1, &dblRealTemp, &dblImgTemp);
295             dblRealTemp                                 = dexps(dblRealTemp * _dblReal2);
296             dblImgTemp                                  = dblImgTemp * _dblReal2;
297             *_pdblRealOut                               = dblRealTemp * dcoss(dblImgTemp);
298             *_pdblImgOut                                = dblRealTemp * dsins(dblImgTemp);
299         }
300         else
301         {
302             if (_dblReal2 > 0)
303             {
304                 //0 ^ R*+
305                 *_pdblRealOut           = 0;
306                 *_pdblImgOut            = 0;
307             }
308             else if (_dblReal2 < 0)
309             {
310                 //0 ^ R*-
311                 //FIXME : ieee
312                 //generate +Inf
313                 double dblZero  = 0.0;
314                 *_pdblRealOut           = 1.0 / (dblZero);
315                 *_pdblImgOut            = 0;
316             }
317             else
318             {
319                 //0 ^ 0
320                 *_pdblRealOut           = 1;
321                 *_pdblImgOut            = 0;
322             }
323         }
324     }
325     return 0;
326 }
327
328 /*wwpowe*/
329 int iPowerComplexScalarByComplexScalar(
330     double _dblReal1, double _dblImg1,
331     double _dblReal2, double _dblImg2,
332     double* _pdblRealOut, double* _pdblImgOut)
333 {
334     if (_dblImg2 == 0)
335     {
336         //C ^ R
337         iPowerComplexScalarByRealScalar(
338             _dblReal1, _dblImg1,
339             _dblReal2,
340             _pdblRealOut, _pdblImgOut);
341     }
342     else
343     {
344         //C ^ C
345         if (dabss(_dblReal1) + dabss(_dblImg1) != 0)
346         {
347             // ! 0 ^ C
348             double dblRealTemp = 0;
349             double dblImgTemp  = 0;
350
351             wlog(_dblReal1, _dblImg1, &dblRealTemp, &dblImgTemp);
352             C2F(wmul)(&dblRealTemp, &dblImgTemp, &_dblReal2, &_dblImg2, &dblRealTemp, &dblImgTemp);
353             dblRealTemp   = dexps(dblRealTemp);
354             *_pdblRealOut = dblRealTemp * dcoss(dblImgTemp);
355             *_pdblImgOut  = dblRealTemp * dsins(dblImgTemp);
356         }
357         else
358         {
359             // 0 ^ C
360             //FIXME : ieee
361             //generate +Inf
362             double dblZero = 0.0;
363             *_pdblRealOut  = 1.0 / (dblZero);
364             *_pdblImgOut   = 0;
365         }
366     }
367     return 0;
368 }
369
370 /*s .^ []*/
371 int iPowerRealScalarByRealMatrix(
372     double _dblReal1,
373     double* _pdblReal2, int _iRows2, int _iCols2,
374     double* _pdblRealOut,       double* _pdblImgOut, int *_piComplex)
375 {
376     int i = 0;
377     *_piComplex = 0;
378     for (i = 0 ; i < _iRows2 * _iCols2 ; i++)
379     {
380         int iComplex = 0;
381         iPowerRealScalarByRealScalar(
382             _dblReal1, _pdblReal2[i], &_pdblRealOut[i], &_pdblImgOut[i], &iComplex);
383         *_piComplex |= iComplex;
384     }
385     return 0;
386 }
387
388
389 int iPowerComplexScalarByRealMatrix(
390     double _dblReal1, double _dblImg1,
391     double* _pdblReal2, int _iRows2, int _iCols2,
392     double* _pdblRealOut,       double* _pdblImgOut)
393 {
394     int i = 0;
395     for (i = 0 ; i < _iRows2 * _iCols2 ; i++)
396     {
397         iPowerComplexScalarByRealScalar(
398             _dblReal1, _dblImg1, _pdblReal2[i], &_pdblRealOut[i], &_pdblImgOut[i]);
399     }
400     return 0;
401 }
402
403 int iPowerRealScalarByComplexMatrix(
404     double _dblReal1,
405     double* _pdblReal2, double* _pdblImg2, int _iRows2, int _iCols2,
406     double* _pdblRealOut,       double* _pdblImgOut)
407 {
408     int i = 0;
409     for (i = 0 ; i < _iRows2 * _iCols2 ; i++)
410     {
411         iPowerRealScalarByComplexScalar(
412             _dblReal1, _pdblReal2[i], _pdblImg2[i], &_pdblRealOut[i], &_pdblImgOut[i]);
413     }
414     return 0;
415 }
416
417 int iPowerComplexScalarByComplexMatrix(
418     double _dblReal1, double _dblImg1,
419     double* _pdblReal2, double* _pdblImg2, int _iRows2, int _iCols2,
420     double* _pdblRealOut,       double* _pdblImgOut)
421 {
422     int i = 0;
423     for (i = 0 ; i < _iRows2 * _iCols2 ; i++)
424     {
425         iPowerComplexScalarByComplexScalar(
426             _dblReal1, _dblImg1, _pdblReal2[i], _pdblImg2[i], &_pdblRealOut[i], &_pdblImgOut[i]);
427     }
428     return 0;
429 }
430
431 //Square Matrix ^ Scalar
432 int iPowerRealSquareMatrixByRealScalar(
433     double* _pdblReal1, int _iRows1, int _iCols1,
434     double _dblReal2,
435     double* _pdblRealOut,       double* _pdblImgOut, int *_iComplex)
436 {
437     int iInv = 0;
438     int iExpRef = (int)_dblReal2;
439     if (iExpRef < 0)
440     {
441         //call matrix invetion
442         iInv = 1;
443         iExpRef = -iExpRef;
444     }
445
446     if ((int)_dblReal2 == _dblReal2) //integer exponent
447     {
448         if (iExpRef == 1)
449         {
450             int iSize = _iRows1 * _iCols1;
451             int iOne = 1;
452             C2F(dcopy)(&iSize, _pdblReal1, &iOne, _pdblRealOut, &iOne);
453         }
454         else if (iExpRef == 0)
455         {
456             int iSize       = _iRows1 * _iCols1;
457             int iOne        = 1;
458             double dblOne   = 1;
459             double dblZero  = 0;
460             int iRowp1      = _iRows1 + 1;
461
462             if (C2F(dasum)(&iSize, _pdblReal1, &iOne) == 0)
463             {
464                 //Invalid exponent
465                 return 1;
466             }
467             C2F(dset)(&iSize, &dblZero, _pdblRealOut, &iOne);
468             C2F(dset)(&_iRows1, &dblOne, _pdblRealOut, &iRowp1);
469         }
470         else
471         {
472             int i           = 0;
473             int iOne        = 1;
474             int iPos        = 0;
475             int iPrevPos    = 0;
476             int iSize       = _iCols1 * _iCols1;
477             int iInc        = _iCols1 + 1;
478
479             double alpha = 1.;
480             double beta  = 0.;
481
482             double* pdblResult  = (double*)malloc(iSize * sizeof(double));
483             double* pdblTmp     = (double*)malloc(iSize * sizeof(double));
484
485             // initialize the output as identity matrix
486             // because we use this array in the first multiplication.
487             memset(_pdblRealOut, 0x00, iSize * sizeof(double));
488             C2F(dset)(&_iCols1, &alpha, _pdblRealOut, &iInc);
489
490             // copy input data to avoid input modification
491             C2F(dcopy)(&iSize, _pdblReal1, &iOne, pdblTmp, &iOne);
492
493             // find all "1" in binary representation of the power integer.
494             // then perform multiplication according the "1" position.
495             while (iExpRef != 0)
496             {
497                 if (iExpRef & 1u)
498                 {
499                     // start the loop at the last position
500                     // because we keep the last result to perform
501                     // the multiplication.
502                     for (i = iPrevPos; i < iPos; i++)
503                     {
504                         double* pdblSwitch;
505                         C2F(dgemm)( "N", "N", &_iCols1, &_iCols1, &_iCols1, &alpha, pdblTmp, &_iCols1,
506                                     pdblTmp, &_iCols1, &beta, pdblResult, &_iCols1);
507
508                         pdblSwitch  = pdblTmp;
509                         pdblTmp     = pdblResult;
510                         pdblResult  = pdblSwitch;
511                     }
512
513                     iPrevPos = iPos;
514
515                     C2F(dgemm)( "N", "N", &_iCols1, &_iCols1, &_iCols1, &alpha, pdblTmp, &_iCols1,
516                                 _pdblRealOut, &_iCols1, &beta, pdblResult, &_iCols1);
517
518                     C2F(dcopy)(&iSize, pdblResult, &iOne, _pdblRealOut, &iOne);
519                 }
520
521                 iPos++;
522                 iExpRef = iExpRef >> 1;
523             }
524
525             free(pdblResult);
526             free(pdblTmp);
527         }
528     }
529     else
530     {
531         //floating point exponent
532         return -1; // manage by overload
533     }
534
535     if (iInv)
536     {
537         double dblRcond;
538         int ret = iInvertMatrixM(_iRows1, _iCols1, _pdblRealOut, 0/* is complex*/, &dblRcond);
539         if (ret == -1)
540         {
541             if (getWarningMode())
542             {
543                 sciprint(_("Warning :\n"));
544                 sciprint(_("matrix is close to singular or badly scaled. rcond = %1.4E\n"), dblRcond);
545                 Sciwarning(_("computing least squares solution. (see lsq).\n"));
546             }
547         }
548     }
549
550     *_iComplex = 0;
551     return 0;
552 }
553
554 int iPowerRealSquareMatrixByComplexScalar(
555     double* _pdblReal1, int _iRows1, int _iCols1,
556     double _dblReal2, double _dblImg2,
557     double* _pdblRealOut,       double* _pdblImgOut)
558 {
559     return 0;
560 }
561
562 int iPowerComplexSquareMatrixByRealScalar(
563     double* _pdblReal1, double* _pdblImg1, int _iRows1, int _iCols1,
564     double _dblReal2,
565     double* _pdblRealOut,       double* _pdblImgOut)
566 {
567     int iInv = 0;
568     int iExpRef = (int)_dblReal2;
569     if (iExpRef < 0)
570     {
571         //call matrix invetion
572         iInv = 1;
573         iExpRef = -iExpRef;
574     }
575
576     if ((int)_dblReal2 == _dblReal2) //integer exponent
577     {
578         if (iExpRef == 1)
579         {
580             int iSize = _iRows1 * _iCols1;
581             int iOne = 1;
582             C2F(dcopy)(&iSize, _pdblReal1, &iOne, _pdblRealOut, &iOne);
583             C2F(dcopy)(&iSize, _pdblImg1, &iOne, _pdblImgOut, &iOne);
584         }
585         else if (iExpRef == 0)
586         {
587             int iSize       = _iRows1 * _iCols1;
588             int iOne        = 1;
589             double dblOne   = 1;
590             double dblZero  = 0;
591             int iRowp1      = _iRows1 + 1;
592
593             if (C2F(dasum)(&iSize, _pdblReal1, &iOne) == 0)
594             {
595                 //Invalid exponent
596                 return 1;
597             }
598             C2F(dset)(&iSize, &dblZero, _pdblRealOut, &iOne);
599             C2F(dset)(&_iRows1, &dblOne, _pdblRealOut, &iRowp1);
600         }
601         else
602         {
603             int i = 0;
604             int iOne = 1;
605             int iTwo = 2;
606             int iPos = 0;
607             int iPrevPos = 0;
608             int iSize = _iCols1 * _iCols1;
609             int iSize2 = iSize * 2;
610             int iInc = (_iCols1 + 1) * 2;
611
612             double alpha[2] = {1., 0.};
613             double beta[2] = {0., 0.};
614
615             double* pdblResult = (double*)malloc(iSize2 * sizeof(double));
616             double* pdblTmp = (double*)malloc(iSize2 * sizeof(double));
617             double* pdblOut = (double*)malloc(iSize2 * sizeof(double));
618
619             // initialize the output as identity matrix
620             // because we use this array in the first multiplication.
621             memset(pdblOut, 0x00, iSize2 * sizeof(double));
622             C2F(dset)(&_iCols1, &alpha[0], pdblOut, &iInc);
623
624             // copy input complex in working array as Z storage.
625             C2F(dcopy)(&iSize, _pdblReal1, &iOne, pdblTmp, &iTwo);
626             C2F(dcopy)(&iSize, _pdblImg1, &iOne, pdblTmp + 1, &iTwo);
627
628             // find all "1" in binary representation of the power integer.
629             // then perform multiplication according the "1" position.
630             while (iExpRef != 0)
631             {
632                 if (iExpRef & 1u)
633                 {
634                     // start the loop at the last position
635                     // because we keep the last result to perform
636                     // the multiplication.
637                     for (i = iPrevPos; i < iPos; i++)
638                     {
639                         double* pdblSwitch;
640                         C2F(zgemm)("N", "N", &_iCols1, &_iCols1, &_iCols1, alpha, pdblTmp, &_iCols1,
641                                    pdblTmp, &_iCols1, beta, pdblResult, &_iCols1);
642
643                         pdblSwitch = pdblTmp;
644                         pdblTmp = pdblResult;
645                         pdblResult = pdblSwitch;
646                     }
647
648                     iPrevPos = iPos;
649
650                     C2F(zgemm)("N", "N", &_iCols1, &_iCols1, &_iCols1, alpha, pdblTmp, &_iCols1,
651                                pdblOut, &_iCols1, beta, pdblResult, &_iCols1);
652
653                     C2F(dcopy)(&iSize2, pdblResult, &iOne, pdblOut, &iOne);
654                 }
655
656                 iPos++;
657                 iExpRef = iExpRef >> 1;
658             }
659
660             C2F(dcopy)(&iSize, pdblOut, &iTwo, _pdblRealOut, &iOne);
661             C2F(dcopy)(&iSize, pdblOut + 1, &iTwo, _pdblImgOut, &iOne);
662
663             free(pdblResult);
664             free(pdblTmp);
665             free(pdblOut);
666         }
667     }
668     else
669     {
670         //floating point exponent
671         return -1; // manage by overload
672     }
673
674     if (iInv)
675     {
676         double dblRcond;
677         double* pData = (double*)oGetDoubleComplexFromPointer(_pdblRealOut, _pdblImgOut, _iRows1 * _iCols1);
678         int ret = iInvertMatrixM(_iRows1, _iCols1, pData, 1/* is complex*/, &dblRcond);
679         if (ret == -1)
680         {
681             if (getWarningMode())
682             {
683                 sciprint(_("Warning :\n"));
684                 sciprint(_("matrix is close to singular or badly scaled. rcond = %1.4E\n"), dblRcond);
685                 Sciwarning(_("computing least squares solution. (see lsq).\n"));
686             }
687         }
688
689         vGetPointerFromDoubleComplex((doublecomplex*)pData, _iRows1 * _iCols1, _pdblRealOut, _pdblImgOut);
690         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
691     }
692
693     return 0;
694 }
695
696 int iPowerComplexSquareMatrixByComplexScalar(
697     double* _pdblReal1, double* _pdblImg1, int _iRows1, int _iCols1,
698     double _dblReal2, double _dblImg2,
699     double* _pdblRealOut,       double* _pdblImgOut)
700 {
701     return 0;
702 }