* Bug 15330 fixed: spec.tst was crashing on Linux
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_spec.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2009 - DIGITEO - Bernard HUGUENEY
4 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14 *
15 */
16 /*--------------------------------------------------------------------------*/
17
18 #include "linear_algebra_gw.hxx"
19 #include "function.hxx"
20 #include "double.hxx"
21 #include "overload.hxx"
22
23 extern "C"
24 {
25 #include "localization.h"
26 #include "Scierror.h"
27 #include "sciprint.h"
28 #include "eigen.h"
29 #include "issymmetric.h"
30 #include "vfinite.h"
31 }
32
33 bool isNoZeroImag(types::Double* _pDbl);
34 /*--------------------------------------------------------------------------*/
35
36 types::Function::ReturnValue sci_spec(types::typed_list &in, int _iRetCount, types::typed_list &out)
37 {
38     double* pDataA          = NULL;
39     double* pDataB          = NULL;
40     bool symmetric          = FALSE;
41     int iRet                = 0;
42
43     _iRetCount = std::max(1, _iRetCount);
44
45     if (in.size() != 1 && in.size() != 2)
46     {
47         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "spec", 1, 2);
48         return types::Function::Error;
49     }
50
51     if (_iRetCount > 2 * in.size())
52     {
53         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "spec", 1, 2 * in.size());
54         return types::Function::Error;
55     }
56
57     if (in[0]->isDouble() == false)
58     {
59         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_spec";
60         return Overload::call(wstFuncName, in, _iRetCount, out);
61     }
62
63     types::Double* in0 = in[0]->getAs<types::Double>();
64
65     if (in0->getCols() != in0->getRows())
66     {
67         Scierror(20, _("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "spec", 1);
68         return types::Function::Error;
69     }
70
71     if (in0->getRows() == -1 || in0->getCols() == -1) // manage eye case
72     {
73         Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "spec", 1);
74         return types::Function::Error;
75     }
76
77     if (in0->getCols() == 0 || in0->getRows() == 0) // size null
78     {
79         out.push_back(types::Double::Empty());
80         for (int i = 1; i < _iRetCount; i++)
81         {
82             out.push_back(types::Double::Empty());
83         }
84         return types::Function::OK;
85     }
86
87     types::Double* pDblA = in0->clone()->getAs<types::Double>();
88
89     if (in.size() == 1)
90     {
91         types::Double* pDblEigenValues  = NULL;
92         types::Double* pDblEigenVectors = NULL;
93
94         if (pDblA->isComplex())
95         {
96             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
97             if (!pDataA)
98             {
99                 pDblA->killMe();
100                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
101                 return types::Function::Error;
102             }
103         }
104         else
105         {
106             pDataA = pDblA->getReal();
107         }
108
109         int totalSize = pDblA->getSize();
110         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
111         {
112             if (pDblA->isComplex())
113             {
114                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
115             }
116             pDblA->killMe();
117             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
118             return types::Function::Error;
119         }
120
121         symmetric = isSymmetric(pDblA->getReal(), pDblA->getImg(), pDblA->isComplex(), pDblA->getRows(), pDblA->getCols()) == 1;
122         int eigenValuesCols = (_iRetCount <= 1) ? 1 : pDblA->getCols();
123
124         if (symmetric)
125         {
126             pDblEigenValues = new types::Double(pDblA->getCols(), eigenValuesCols);
127             if (_iRetCount == 2)
128             {
129                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), pDblA->isComplex());
130             }
131         }
132         else
133         {
134             pDblEigenValues  = new types::Double(pDblA->getCols(), eigenValuesCols, true);
135             if (_iRetCount == 2)
136             {
137                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), true);
138             }
139         }
140
141         if (pDblA->isComplex())
142         {
143             if (symmetric)
144             {
145                 iRet = iEigen1ComplexSymmetricM((doublecomplex*)pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
146
147                 if (iRet < 0)
148                 {
149                     vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
150                     pDblA->killMe();
151                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec");
152                     return types::Function::Error;
153                 }
154
155                 if (iRet > 0)
156                 {
157                     vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
158                     pDblA->killMe();
159                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
160                     return types::Function::Error;
161                 }
162
163                 if (_iRetCount == 2)
164                 {
165                     vGetPointerFromDoubleComplex((doublecomplex*)pDataA, pDblA->getSize(), pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
166                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
167                     out.push_back(pDblEigenVectors);
168                 }
169                 out.push_back(pDblEigenValues);
170                 pDblA->killMe();
171                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
172             }
173             else // not symmetric
174             {
175                 doublecomplex* pEigenValues = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
176                 doublecomplex* pEigenVectors = pDblEigenVectors ? (doublecomplex*)MALLOC(sizeof(doublecomplex) * pDblA->getSize()) : NULL;
177                 iRet = iEigen1ComplexM((doublecomplex*)pDataA, pDblA->getCols(), pEigenValues, pEigenVectors);
178                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
179                 if (iRet < 0)
180                 {
181                     pDblA->killMe();
182                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec");
183                     return types::Function::Error;
184                 }
185
186                 if (iRet > 0)
187                 {
188                     pDblA->killMe();
189                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
190                     return types::Function::Error;
191                 }
192
193                 if (_iRetCount == 2)
194                 {
195                     expandZToDiagonalOfCMatrix(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
196                     vGetPointerFromDoubleComplex(pEigenVectors, pDblA->getSize(), pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
197
198                     FREE(pEigenVectors);
199                     out.push_back(pDblEigenVectors);
200                 }
201                 else
202                 {
203                     vGetPointerFromDoubleComplex(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
204                 }
205                 out.push_back(pDblEigenValues);
206                 FREE(pEigenValues);
207                 pDblA->killMe();
208             }
209         }
210         else // real
211         {
212             if (symmetric)
213             {
214                 iRet = iEigen1RealSymmetricM(pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
215                 if (iRet < 0)
216                 {
217                     pDblA->killMe();
218                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec");
219                     return types::Function::Error;
220                 }
221
222                 if (iRet > 0)
223                 {
224                     pDblA->killMe();
225                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
226                     return types::Function::Error;
227                 }
228
229                 if (_iRetCount == 2)
230                 {
231                     if (pDblEigenVectors)
232                     {
233                         pDblEigenVectors->killMe();
234                     }
235                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
236                     out.push_back(pDblA);
237                 }
238                 else
239                 {
240                     pDblA->killMe();
241                 }
242
243                 out.push_back(pDblEigenValues);
244             }
245             else // not symmetric
246             {
247                 iRet = iEigen1RealM(pDataA, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg(), pDblEigenVectors ? pDblEigenVectors->getReal() : NULL, pDblEigenVectors ? pDblEigenVectors->getImg() : NULL);
248
249                 if (iRet < 0)
250                 {
251                     pDblA->killMe();
252                     if (pDblEigenVectors)
253                     {
254                         pDblEigenVectors->killMe();
255                     }
256                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec");
257                     return types::Function::Error;
258                 }
259
260                 if (iRet > 0)
261                 {
262                     pDblA->killMe();
263                     if (pDblEigenVectors)
264                     {
265                         pDblEigenVectors->killMe();
266                     }
267                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of WR and WI contain eigenvalues which have converged.\n"), "spec", iRet);
268                     return types::Function::Error;
269                 }
270
271                 if (_iRetCount == 2)
272                 {
273                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
274                     expandToDiagonalOfMatrix(pDblEigenValues->getImg(), pDblA->getCols());
275                     out.push_back(pDblEigenVectors);
276                 }
277
278                 out.push_back(pDblEigenValues);
279                 pDblA->killMe();
280             }
281         }
282
283         return types::Function::OK;
284     }
285
286     if (in.size() == 2)
287     {
288         types::Double* pDblL            = NULL;
289         types::Double* pDblR            = NULL;
290         types::Double* pDblBeta         = NULL;
291         types::Double* pDblAlpha        = NULL;
292         doublecomplex* pL               = NULL;
293         doublecomplex* pR               = NULL;
294         doublecomplex* pBeta            = NULL;
295         doublecomplex* pAlpha           = NULL;
296         bool bIsComplex                 = false;
297
298         if (in[1]->isDouble() == false)
299         {
300             std::wstring wstFuncName = L"%" + in[1]->getShortTypeStr() + L"_spec";
301             return Overload::call(wstFuncName, in, _iRetCount, out);
302         }
303
304         types::Double* in1 = in[1]->getAs<types::Double>();
305
306         if (in1->getCols() != in1->getRows())
307         {
308             Scierror(20, _("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "spec", 2);
309             return types::Function::Error;
310         }
311
312         if (pDblA->getRows() != in1->getRows() && pDblA->getCols() != in1->getCols())
313         {
314             pDblA->killMe();
315             Scierror(999, _("%s: Arguments %d and %d must have equal dimensions.\n"), "spec", 1, 2);
316             return types::Function::Error;
317         }
318
319         //chekc if A and B are real complex or with imag part at 0
320         if (isNoZeroImag(pDblA) == false && isNoZeroImag(in1) == false)
321         {
322             //view A and B as real matrix
323             bIsComplex = false;
324         }
325         else
326         {
327             bIsComplex = pDblA->isComplex() || in1->isComplex();
328         }
329
330         types::Double* pDblB = in1->clone()->getAs<types::Double>();
331         if (bIsComplex)
332         {
333             if (pDblA->isComplex() == false)
334             {
335                 pDblA->setComplex(true);
336             }
337
338             if (pDblB->isComplex() == false)
339             {
340                 pDblB->setComplex(true);
341             }
342
343             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
344             pDataB = (double*)oGetDoubleComplexFromPointer(pDblB->getReal(), pDblB->getImg(), pDblB->getSize());
345
346             if (!pDataA && !pDataB)
347             {
348                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
349                 return types::Function::Error;
350             }
351
352             if (!pDataA)
353             {
354                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
355                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
356                 return types::Function::Error;
357             }
358
359             if (!pDataB)
360             {
361                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
362                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
363                 return types::Function::Error;
364             }
365         }
366         else
367         {
368             pDataA = pDblA->getReal();
369             pDataB = pDblB->getReal();
370         }
371
372         int totalSize = pDblA->getSize();
373
374         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
375         {
376             pDblA->killMe();
377             pDblB->killMe();
378             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
379             return types::Function::Error;
380         }
381
382         if ((pDblB->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataB) : C2F(vfinite)(&totalSize, pDataB)) == false)
383         {
384             pDblA->killMe();
385             pDblB->killMe();
386             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 2);
387             return types::Function::Error;
388         }
389
390         switch (_iRetCount)
391         {
392             case 4:
393             {
394                 pDblL = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
395                 if (bIsComplex)
396                 {
397                     pL = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
398                 }
399             }
400             case 3:
401             {
402                 pDblR = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
403                 if (bIsComplex)
404                 {
405                     pR = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
406                 }
407             }
408             case 2:
409             {
410                 if (bIsComplex)
411                 {
412                     pBeta = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
413                 }
414                 pDblBeta = new types::Double(pDblA->getCols(), 1, pBeta ? true : false);
415             }
416             default : // case 1:
417             {
418                 if (bIsComplex)
419                 {
420                     pAlpha = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
421                 }
422                 pDblAlpha = new types::Double(pDblA->getCols(), 1, true);
423             }
424         }
425
426         if (bIsComplex)
427         {
428             iRet = iEigen2ComplexM((doublecomplex*)pDataA, (doublecomplex*)pDataB, pDblA->getCols(), pAlpha, pBeta, pR, pL);
429         }
430         else
431         {
432             iRet = iEigen2RealM(    pDataA, pDataB, pDblA->getCols(),
433                                     pDblAlpha->getReal(), pDblAlpha->getImg(),
434                                     pDblBeta ? pDblBeta->getReal()  : NULL,
435                                     pDblR    ? pDblR->getReal()     : NULL,
436                                     pDblR    ? pDblR->getImg()      : NULL,
437                                     pDblL    ? pDblL->getReal()     : NULL,
438                                     pDblL    ? pDblL->getImg()      : NULL);
439         }
440
441         if (iRet > 0)
442         {
443             sciprint(_("Warning :\n"));
444             sciprint(_("Non convergence in the QZ algorithm.\n"));
445             sciprint(_("The top %d  x %d blocks may not be in generalized Schur form.\n"), iRet);
446         }
447
448         if (iRet < 0)
449         {
450             pDblA->killMe();
451             pDblB->killMe();
452             Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec");
453             return types::Function::Error;
454         }
455
456         if (iRet > 0)
457         {
458             if (bIsComplex)
459             {
460                 if (iRet <= pDblA->getCols())
461                 {
462                     Scierror(24, _("%s: The QZ iteration failed in DGGEV.\n"), "spec");
463                 }
464                 else
465                 {
466                     if (iRet == pDblA->getCols() + 1)
467                     {
468                         Scierror(999, _("%s: Other than QZ iteration failed in DHGEQZ.\n"), "spec");
469                     }
470                     if (iRet == pDblA->getCols() + 2)
471                     {
472                         Scierror(999, _("%s: Error return from DTGEVC.\n"), "spec");
473                     }
474                 }
475             }
476             else
477             {
478                 Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
479             }
480
481             pDblA->killMe();
482             pDblB->killMe();
483             if (pDataA)
484             {
485                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
486             }
487
488             if (pDataB)
489             {
490                 vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
491             }
492             return types::Function::Error;
493         }
494
495         if (bIsComplex)
496         {
497             switch (_iRetCount)
498             {
499                 case 4:
500                     vGetPointerFromDoubleComplex(pL, pDblA->getSize(), pDblL->getReal(), pDblL->getImg());
501                 case 3:
502                     vGetPointerFromDoubleComplex(pR, pDblA->getSize(), pDblR->getReal(), pDblR->getImg());
503                 case 2:
504                     vGetPointerFromDoubleComplex(pBeta, pDblA->getCols(), pDblBeta->getReal(), pDblBeta->getImg());
505                 default : // case 1:
506                     vGetPointerFromDoubleComplex(pAlpha, pDblA->getCols(), pDblAlpha->getReal(), pDblAlpha->getImg());
507             }
508         }
509
510         switch (_iRetCount)
511         {
512             case 1:
513             {
514                 out.push_back(pDblAlpha);
515                 break;
516             }
517             case 2:
518             {
519                 out.push_back(pDblAlpha);
520                 out.push_back(pDblBeta);
521                 break;
522             }
523             case 3:
524             {
525                 out.push_back(pDblAlpha);
526                 out.push_back(pDblBeta);
527                 out.push_back(pDblR);
528                 break;
529             }
530             case 4:
531             {
532                 out.push_back(pDblAlpha);
533                 out.push_back(pDblBeta);
534                 out.push_back(pDblL);
535                 out.push_back(pDblR);
536             }
537         }
538
539         if (pAlpha)
540         {
541             vFreeDoubleComplexFromPointer(pAlpha);
542         }
543         if (pBeta)
544         {
545             vFreeDoubleComplexFromPointer(pBeta);
546         }
547         if (pL)
548         {
549             vFreeDoubleComplexFromPointer(pL);
550         }
551         if (pR)
552         {
553             vFreeDoubleComplexFromPointer(pR);
554         }
555         if (bIsComplex && pDblB->isComplex())
556         {
557             vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
558         }
559         pDblB->killMe();
560
561     } // if(in.size() == 2)
562
563     if (pDblA->isComplex())
564     {
565         vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
566     }
567
568     return types::Function::OK;
569 }
570 /*--------------------------------------------------------------------------*/
571 bool isNoZeroImag(types::Double* _pDbl)
572 {
573     double* pdbl = _pDbl->getImg();
574     if (pdbl)
575     {
576         for (int i = 0 ; i < _pDbl->getSize() ; i++)
577         {
578             if (pdbl[i])
579             {
580                 return true;
581             }
582         }
583     }
584     return false;
585 }
586 /*--------------------------------------------------------------------------*/