linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_chol.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 "chol.h"
26 }
27 /*--------------------------------------------------------------------------*/
28
29 types::Function::ReturnValue sci_chol(types::typed_list &in, int _iRetCount, types::typed_list &out)
30 {
31     types::Double* pDbl     = NULL;
32     double* pdOut           = NULL;
33     int iCholProductResult  = 0;
34     types::Double* result   = NULL;
35
36     if(in.size() != 1)
37     {
38         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"chol", 1);
39         return types::Function::Error;
40     }
41
42     if((in[0]->isDouble() == false))
43     {
44         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_chol";
45         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
46     }
47
48     pDbl = in[0]->getAs<types::Double>();
49     if(pDbl->getRows() != pDbl->getCols())
50     {
51                 ScierrorW(20, _W("%ls: Wrong type for argument %d: Square matrix expected.\n"), L"chol", 1);
52         return types::Function::Error;            
53     }
54
55     if(pDbl->getRows() == 0)
56     {
57         out.push_back(types::Double::Empty());
58         return types::Function::OK;
59     }
60     else if(pDbl->getRows() == -1) // manage eye case
61     {
62         if(pDbl->get(0) <= 0)
63         {
64                 /* Matrix must be positive definite */
65                         ScierrorW(29, _W("%ls: Matrix is not positive definite.\n"), L"chol");
66             return types::Function::Error;
67         }
68
69         types::Double* result = new types::Double(sqrt(pDbl->get(0)));
70         out.push_back(result);
71         return types::Function::OK;
72     }
73
74     if(pDbl->isComplex())
75     {
76         doublecomplex* poData = oGetDoubleComplexFromPointer(pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
77                 iCholProductResult = iComplexCholProduct(poData, pDbl->getRows());
78
79         result = new types::Double(pDbl->getRows(), pDbl->getCols(), true);
80         
81                 vGetPointerFromDoubleComplex(poData, pDbl->getSize(), result->getReal(), result->getImg());
82                 vFreeDoubleComplexFromPointer(poData);
83     }
84     else
85     {
86         result = pDbl->clone()->getAs<types::Double>();
87                 iCholProductResult = iRealCholProduct(result->get(), pDbl->getRows());
88     }
89
90     if(iCholProductResult > 0)
91     {
92         /* Matrix must be positive definite */
93                 ScierrorW(29, _W("%ls: Matrix is not positive definite.\n"), L"chol");
94         return types::Function::Error;
95     }
96
97     out.push_back(result);
98     return types::Function::OK;
99 }
100 /*--------------------------------------------------------------------------*/
101