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