linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_inv.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 "invert_matrix.h"
27 #include "sci_warning.h"
28 }
29 /*--------------------------------------------------------------------------*/
30
31 types::Function::ReturnValue sci_inv(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDbl = NULL;
34     double* pData       = NULL;
35     int ret             = 0;
36
37     if(in.size() != 1)
38     {
39         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"inv", 1);
40         return types::Function::Error;
41     }
42
43     if((in[0]->isDouble() == false))
44     {
45         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_inv";
46         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
47     }
48
49     pDbl = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>(); // input data will be modified
50
51     if(pDbl->getRows() != pDbl->getCols())
52     {
53         ScierrorW(20, _W("%ls: Wrong type for argument %d: Square matrix expected.\n"), L"inv", 1);
54         return types::Function::Error;      
55     }
56
57     if(pDbl->getRows() == 0)
58     {
59         out.push_back(types::Double::Empty());
60         return types::Function::OK;
61     }
62
63     if(pDbl->isComplex())
64     {
65         /* c -> z */
66         pData = (double*)oGetDoubleComplexFromPointer( pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
67     }
68     else
69     {
70         pData = pDbl->getReal();
71     }
72
73     if(pDbl->getCols() == -1)
74     {
75         pData[0] = 1./pData[0];
76     }
77     else
78     {
79         double dblRcond;
80         ret = iInvertMatrixM(pDbl->getRows(), pDbl->getCols(), pData, pDbl->isComplex(), &dblRcond);
81         if(pDbl->isComplex())
82         {
83             /* z -> c */
84             vGetPointerFromDoubleComplex((doublecomplex*)pData, pDbl->getSize(), pDbl->getReal(), pDbl->getImg());
85             vFreeDoubleComplexFromPointer((doublecomplex*)pData);
86         }
87
88         if(ret == -1)
89         {
90             if(getWarningMode())
91             {
92                 sciprint(_("Warning :\n"));
93                 sciprint(_("matrix is close to singular or badly scaled. rcond = %1.4E\n"), dblRcond);
94                 sciprint(_("computing least squares solution. (see lsq).\n"));
95             }
96             }
97     }
98
99     if(ret == 19)
100     {
101         ScierrorW(19, _W("%ls: Problem is singular.\n"), L"inv");
102         return types::Function::Error;
103     }
104
105     out.push_back(pDbl);
106     return types::Function::OK;
107 }
108 /*--------------------------------------------------------------------------*/
109