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