linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_bdiag.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
19 extern "C"
20 {
21 #include "localization.h"
22 #include "Scierror.h"
23 #include "core_math.h" /* for Abs  macro */
24 }
25 /*--------------------------------------------------------------------------*/
26 extern "C"
27 {
28 extern int C2F(vfinite)(int const*n, double const*v);
29
30 extern void C2F(bdiag)(int const * lda, int const* n, double* a, double const * epsshr, double const* rMax, double* er, double* ei, int *bs
31                        , double* x, double* xi, double* scale, int const* job, int * fail);
32 extern void C2F(wbdiag)(int const* lda, int const* n, double* ar, double* ai, double const*rMax, double* er, double* ei, int* bs
33                         , double* xr , double* xi, double* yr, double* yi, double* scale, int const* job, int* fail);
34 }
35 /*--------------------------------------------------------------------------*/
36 types::Function::ReturnValue sci_bdiag(types::typed_list &in, int _iRetCount, types::typed_list &out)
37 {
38     types::Double* pDblMatrix   = NULL;
39     types::Double* pDblScalar   = NULL;
40     double rMax                 = 0;
41
42     if((in.size() != 1) && (in.size() != 2))
43     {
44         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d to %d expected.\n"), L"bdiag", 1, 2);
45         return types::Function::Error;
46     }
47     if(_iRetCount > 3)
48     {
49         ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"bdiag", 1, 3);
50         return types::Function::Error;
51     }
52
53     if(in[0]->isDouble() == false)
54     {
55         ScierrorW(201, _W("%ls: Wrong type for argument %d: Real or complex matrix expected.\n"), L"bdiag", 1);
56         return types::Function::Error;     
57     }
58
59     pDblMatrix = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>(); // input data will be modified
60     if(pDblMatrix->getRows() != pDblMatrix->getCols())
61     {
62                 ScierrorW(20, _W("%ls: Wrong type for argument %d: Square matrix expected.\n"), L"bdiag", 1);
63         return types::Function::Error;      
64     }
65
66     if(in.size() == 2)
67     {
68         if(in[1]->isDouble() == false && in[1]->getAs<types::Double>()->isScalar() == false)
69         {
70             ScierrorW(999, _W("%ls: Wrong type for argument %d: A scalar expected.\n"), L"bdiag", 2);
71             return types::Function::Error;     
72         }
73
74         pDblScalar = in[1]->getAs<types::Double>();
75     }
76
77     if(pDblMatrix->getCols() == 0)
78     {
79         int value = Min(_iRetCount, 3);
80         for(int i = 0; i < value; i++)
81         {
82             out.push_back(types::Double::Empty());
83         }
84         return types::Function::OK;
85     }
86
87     int const totalSize = pDblMatrix->getSize();
88
89     if( C2F(vfinite)(&totalSize, pDblMatrix->getReal()) == false &&
90         (pDblMatrix->isComplex() == false || 
91          C2F(vfinite)(&totalSize, pDblMatrix->getImg())))
92     {
93                 ScierrorW(264, _W("%ls: Wrong value for argument %d: Must not contain NaN or Inf.\n"), L"bdiag", 1);
94         return types::Function::Error;   
95     }
96
97     if(pDblScalar != NULL)
98     {
99         rMax = pDblScalar->get(0);
100     }
101     else
102     {
103         for(int j = 0; j < pDblMatrix->getCols(); j++)
104         {
105             double t = 0.0;
106             for(int i = 0; i < pDblMatrix->getCols(); i++)
107             {
108                 t += Abs(pDblMatrix->get(i+j*pDblMatrix->getCols()));
109             }
110             rMax = Max(t, rMax);
111         }
112     }
113
114     types::Double* lx = pDblMatrix->clone()->getAs<types::Double>();
115     int iDim        = pDblMatrix->getCols(); // square matrix
116     int fail        = 0;
117     int const job   = 0;
118
119     /* allocating the two memory buffers in one place as the original code did */
120     double* le = (double*) MALLOC( 2*iDim * sizeof(double) );
121     int*    lib= (int*) MALLOC(iDim * sizeof(int));
122     double* lw = (double*)MALLOC(iDim * sizeof(double));
123
124     if((le && lib && lw) == false)
125     {
126                 ScierrorW(999, _W("%ls: Allocation failed.\n"), L"bdiag");
127         return types::Function::Error;   
128     }
129
130     if(pDblMatrix->isComplex())
131     {
132         double dummy;
133         C2F(wbdiag)(&iDim, &iDim, pDblMatrix->getReal(), pDblMatrix->getImg(), &rMax, le, le + iDim, lib, lx->getReal() , lx->getImg(),  &dummy, &dummy, lw, &job, &fail);
134     }
135     else
136     {
137         double const epsshr= 0.;
138         C2F(bdiag)(&iDim, &iDim, pDblMatrix->getReal(), &epsshr, &rMax, le, le + iDim, lib, lx->getReal(), NULL /* not used when job = 0 */, lw, &job, &fail);
139     }
140
141     if(fail)
142     {
143         FREE(le);
144         FREE(lib);
145         FREE(lw);
146                 ScierrorW(24, _W("%ls: Non convergence in QR steps.\n"), L"bdiag");
147         return types::Function::Error;   
148     }
149     else
150     {
151         out.push_back(pDblMatrix);
152         if(_iRetCount > 1)
153         {
154             out.push_back(lx);
155         }
156         if(_iRetCount == 3)
157         {
158             int nbloc   = 0;
159             int k       = 0;
160
161             for(k = 0; k < iDim; k++)
162             {
163                 if(lib[k] >= 0)
164                 {
165                     nbloc++;
166                 }
167             }
168
169             types::Double* lbs = new types::Double(nbloc, 1);
170
171             for (int i = k = 0; k < iDim; k++)
172             {
173                 if(lib[k] >= 0)
174                 {
175                     lbs->set(i,(double) lib[k]);
176                     i++;
177                 }
178             }
179             out.push_back(lbs);
180         }
181         FREE(le);
182         FREE(lib);
183         FREE(lw);
184     }
185     return types::Function::OK;
186 }
187 /*--------------------------------------------------------------------------*/
188