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