scilab/modules/sparse/sci_gateway/cpp/sci_lusolve.cpp

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) 2014 - Scilab Enterprises - Anais AUBERT
   4  *
   5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
   6  *
   7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
   8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
   9  * This file was originally licensed under the terms of the CeCILL v2.1,
  10  * and continues to be available under such terms.
  11  * For more information, see the COPYING file which you should have received
  12  * along with this program.
  13  *
  14  */
  15
  16 #include "sparse_gw.hxx"
  17 #include "function.hxx"
  18 #include "types.hxx"
  19 #include "double.hxx"
  20 #include "sparse.hxx"
  21 #include "pointer.hxx"
  22 #include "overload.hxx"
  23
  24 extern "C"
  25 {
  26 #include "Scierror.h"
  27 #include "localization.h"
  28 #include "elem_common.h"
  29 #include "lu.h"
  30 }
  31
  32 types::Function::ReturnValue sci_lusolve(types::typed_list &in, int _iRetCount, types::typed_list &out)
  33 {
  34     double abstol   = 0;
  35     double reltol   = 0.001;
  36     int nrank       = 0;
  37     int ierr        = 0;
  38     int m1          = 0;
  39     int n1          = 0;
  40     int m2          = 0;
  41     int n2          = 0;
  42     int nonZeros    = 0;
  43     bool fact       = false;
  44
  45     const void *pData   = NULL;
  46     double* dbl         = NULL;
  47     int* colPos         = NULL;
  48     int* itemsRow       = NULL;
  49     int *fmatindex      = NULL;
  50
  51     //check input parameters
  52     if (in.size() > 2)
  53     {
  54         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "lusolve", 2);
  55         return types::Function::Error;
  56     }
  57
  58     if (_iRetCount > 1)
  59     {
  60         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "lusolve", 1);
  61         return types::Function::Error;
  62     }
  63
  64     if (in[0]->isPointer())
  65     {
  66         types::Pointer *pPointerIn = in[0]->getAs<types::Pointer>();
  67         m1 = pPointerIn->getRows();
  68         n1 = pPointerIn->getCols();
  69         pData = pPointerIn->get();
  70         fmatindex = (int*)pData;
  71         fact = false;
  72     }
  73     else if (in[0]->isSparse())
  74     {
  75         types::Sparse *pSpIn = in[0]->getAs<types::Sparse>();
  76         m1 = pSpIn->getRows();
  77         n1 = pSpIn->getCols();
  78         if (m1 != n1)
  79         {
  80             Scierror(77, _("%s: Wrong size for input argument #%d: Square matrix expected.\n"), "lusolve", 1);
  81             return types::Function::Error;
  82         }
  83
  84         if (pSpIn->isComplex())
  85         {
  86             Scierror(77, _("%s: Wrong type for argument %d: Real matrix expected.\n"), "lusolve", 1);
  87             return types::Function::Error;
  88         }
  89
  90         nonZeros = (int)pSpIn->nonZeros();
  91         dbl = new double[nonZeros];
  92         pSpIn->outputValues(dbl, NULL);
  93         colPos = new int[nonZeros];
  94         itemsRow = new int[m1];
  95         pSpIn->getColPos(colPos);
  96         pSpIn->getNbItemByRow(itemsRow);
  97
  98         fmatindex = new int[1];
  99         abstol = nc_eps_machine();
 100         C2F(lufact1)(dbl, itemsRow, colPos, &m1, &nonZeros, fmatindex, &abstol, &reltol, &nrank, &ierr);
 101         fact = true;
 102
 103         delete[] dbl;
 104         delete[] colPos;
 105         delete[] itemsRow;
 106     }
 107     else
 108     {
 109         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lusolve";
 110         return Overload::call(wstFuncName, in, _iRetCount, out);
 111     }
 112
 113     if ((in[1]->isSparse() == false) && (in[1]->isDouble() == false))
 114     {
 115         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lusolve";
 116         return Overload::call(wstFuncName, in, _iRetCount, out);
 117     }
 118
 119     if (in[1]->isSparse() )
 120     {
 121         Scierror(999, _("%s not yet implemented for full input parameter.\n"), "lusolve");
 122         return types::Function::Error;
 123     }
 124
 125     if (in[1]->isDouble() )
 126     {
 127         types::Double *pDblIn = in[1]->getAs<types::Double>();
 128
 129         m2 = pDblIn->getRows();
 130         n2 = pDblIn->getCols();
 131
 132         if (m2 != m1)
 133         {
 134             Scierror(999, _("%s: Wrong size for input argument #%d: Incompatible dimensions.\n"), "lusolve", 2);
 135             return types::Function::Error;
 136         }
 137
 138         double *dbl  = pDblIn->getReal();
 139         types::Double *pDblOut = new types::Double(m2, n2, pDblIn->isComplex());
 140         double *oReal = pDblOut->get();
 141
 142         if (pDblIn->isComplex())
 143         {
 144             double *imag = pDblIn->getImg();
 145             double *oImg = pDblOut->getImg();
 146             for (int i = 0; i < n2; i++)
 147             {
 148                 int iPos = i * m2;
 149                 C2F(lusolve1)(fmatindex, dbl + iPos, oReal + iPos, &ierr);
 150                 if (ierr > 0)
 151                 {
 152                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
 153                     return types::Function::Error;
 154                 }
 155                 C2F(lusolve1)(fmatindex, imag + iPos, oImg + iPos, &ierr);
 156                 if (ierr > 0)
 157                 {
 158                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
 159                     return types::Function::Error;
 160                 }
 161             }
 162         }
 163         else
 164         {
 165             for (int i = 0; i < n2; i++)
 166             {
 167                 C2F(lusolve1)(fmatindex, &dbl[i * m2], &oReal[i * m2], &ierr);
 168                 if (ierr > 0)
 169                 {
 170                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
 171                     return types::Function::Error;
 172                 }
 173             }
 174         }
 175
 176         if (fact)
 177         {
 178             C2F(ludel1)(fmatindex, &ierr);
 179         }
 180
 181         out.push_back(pDblOut);
 182     }
 183
 184     return types::Function::OK;
 185 }