scilab/modules/cacsd/sci_gateway/cpp/sci_residu.cpp

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) 2014 - Scilab Enterprises - Cedric DELAMARRE
   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 "cacsd_gw.hxx"
  17 #include "function.hxx"
  18 #include "double.hxx"
  19 #include "polynom.hxx"
  20
  21 extern "C"
  22 {
  23 #include "Scierror.h"
  24 #include "localization.h"
  25 #include "elem_common.h"
  26
  27     extern void C2F(residu)(double*, int*, double*, int*, double*, int*, double*, double*, int*);
  28     extern void C2F(wesidu)(double*, double*, int*, double*, double*, int*,
  29                             double*, double*, int*, double*, double*, double*, int*);
  30 }
  31
  32 /*--------------------------------------------------------------------------*/
  33 types::Function::ReturnValue sci_residu(types::typed_list &in, int _iRetCount, types::typed_list &out)
  34 {
  35     int iRows[3]             = {0, 0, 0};
  36     int iCols[3]             = {0, 0, 0};
  37     int iComplex[3]          = {0, 0, 0};
  38     int* piRank[3]           = {NULL, NULL, NULL};
  39     double** pdblInR[3]      = {NULL, NULL, NULL};
  40     double** pdblInI[3]      = {NULL, NULL, NULL};
  41     bool isDeletable[3]      = {false, false, false};
  42     types::Double* pDblIn[3] = {NULL, NULL, NULL};
  43     types::Polynom* pPoly[3] = {NULL, NULL, NULL};
  44     types::Double* pDblOut   = NULL;
  45
  46     double dblEps   = nc_eps_machine();
  47     double dZero    = 0;
  48     int iOne        = 1;
  49     int iSize       = 0;
  50     int iError      = 0;
  51
  52     if (in.size() != 3)
  53     {
  54         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "residu", 3);
  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"), "residu", 1);
  61         return types::Function::Error;
  62     }
  63
  64     try
  65     {
  66         /*** get inputs arguments ***/
  67         for (int i = 0; i < 3; i++)
  68         {
  69             if (in[i]->isDouble())
  70             {
  71                 pDblIn[i] = in[i]->clone()->getAs<types::Double>();
  72                 iRows[i] = pDblIn[i]->getRows();
  73                 iCols[i] = pDblIn[i]->getCols();
  74
  75                 iSize = pDblIn[i]->getSize();
  76                 piRank[i] = new int[iSize];
  77                 memset(piRank[i], 0x00, iSize * sizeof(int));
  78
  79                 pdblInR[i]  = new double*[iSize];
  80                 double* pdbl = pDblIn[i]->get();
  81                 for (int j = 0; j < iSize; j++)
  82                 {
  83                     pdblInR[i][j] = pdbl + j;
  84                 }
  85
  86                 if (pDblIn[i]->isComplex())
  87                 {
  88                     pdblInI[i]  = new double*[iSize];
  89                     double* pdbl = pDblIn[i]->get();
  90                     for (int j = 0; j < iSize; j++)
  91                     {
  92                         pdblInI[i][j] = pdbl + j;
  93                     }
  94                 }
  95             }
  96             else if (in[i]->isPoly())
  97             {
  98                 pPoly[i] = in[i]->clone()->getAs<types::Polynom>();
  99                 iRows[i] = pPoly[i]->getRows();
 100                 iCols[i] = pPoly[i]->getCols();
 101
 102                 iSize = pPoly[i]->getSize();
 103                 piRank[i] = new int[iSize];
 104                 pPoly[i]->getRank(piRank[i]);
 105
 106                 pdblInR[i] = new double*[iSize];
 107                 if (pPoly[i]->isComplex())
 108                 {
 109                     pdblInI[i] = new double*[iSize];
 110                     for (int j = 0; j < iSize; j++)
 111                     {
 112                         pdblInR[i][j] = pPoly[i]->get(j)->get();
 113                         pdblInI[i][j] = pPoly[i]->get(j)->getImg();
 114                     }
 115                 }
 116                 else
 117                 {
 118                     for (int j = 0; j < iSize; j++)
 119                     {
 120                         pdblInR[i][j] = pPoly[i]->get(j)->get();
 121                     }
 122                 }
 123             }
 124             else
 125             {
 126                 Scierror(999, _("%s: Wrong type for input argument #%d: A Matrix or polynom expected.\n"), "residu", i + 1);
 127                 throw 1;
 128             }
 129         }
 130
 131         if (iRows[0] != iRows[1] || iCols[0] != iCols[1] || iRows[0] != iRows[2] || iCols[0] != iCols[2])
 132         {
 133             Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "residu");
 134             throw 1;
 135         }
 136
 137         /*** perform operations ***/
 138         if (pdblInI[0] == NULL && pdblInI[1] == NULL && pdblInI[2] == NULL)
 139         {
 140             // real case
 141             pDblOut = new types::Double(iRows[0], iCols[0]);
 142             double* pdblOut = pDblOut->get();
 143             for (int i = 0; i < iRows[0] * iCols[0]; i++)
 144             {
 145                 int iErr = 0;
 146                 double v = 0;
 147
 148                 C2F(residu)(pdblInR[0][i], piRank[0]+i, pdblInR[1][i], piRank[1]+i,
 149                             pdblInR[2][i], piRank[2]+i, &v, &dblEps, &iErr);
 150                 if (iErr)
 151                 {
 152                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "residu");
 153                     throw iErr;
 154                 }
 155
 156                 pdblOut[i] = v;
 157             }
 158         }
 159         else
 160         {
 161             // complex case
 162             pDblOut = new types::Double(iRows[0], iCols[0], true);
 163             double* pdblOutR = pDblOut->get();
 164             double* pdblOutI = pDblOut->getImg();
 165
 166             for (int i = 0; i < 3; i++)
 167             {
 168                 if (pdblInI[i] == NULL)
 169                 {
 170                     pdblInI[i] = new double*[iSize];
 171                     for (int j = 0; j < iSize; j++)
 172                     {
 173                         int iLen = piRank[i][j] + 1;
 174                         pdblInI[i][j] = new double[iLen];
 175                         memset(pdblInI[i][j], 0x00, iLen * sizeof(double));
 176                     }
 177
 178                     isDeletable[i] = true;
 179                 }
 180             }
 181
 182             for (int i = 0; i < iRows[0] * iCols[0]; i++)
 183             {
 184                 int iErr    = 0;
 185                 double real = 0;
 186                 double imag = 0;
 187
 188                 C2F(wesidu)(pdblInR[0][i], pdblInI[0][i], (piRank[0]) + i,
 189                             pdblInR[1][i], pdblInI[1][i], (piRank[1]) + i,
 190                             pdblInR[2][i], pdblInI[2][i], (piRank[2]) + i,
 191                             &real, &imag, &dblEps, &iErr);
 192
 193                 if (iErr)
 194                 {
 195                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "wesidu");
 196                     throw iErr;
 197                 }
 198
 199                 pdblOutR[i] = real;
 200                 pdblOutI[i] = imag;
 201             }
 202         }
 203     }
 204     catch (int error)
 205     {
 206         iError = error;
 207     }
 208
 209     // free memory
 210     for (int i = 0; i < 3; i++)
 211     {
 212         if (pDblIn[i])
 213         {
 214             delete pDblIn[i];
 215         }
 216
 217         if (pPoly[i])
 218         {
 219             delete pPoly[i];
 220         }
 221
 222         if (piRank[i])
 223         {
 224             delete[] piRank[i];
 225         }
 226
 227         if (pdblInR[i])
 228         {
 229             delete[] pdblInR[i];
 230         }
 231
 232         if (isDeletable[i])
 233         {
 234             for (int j = 0; j < iSize; j++)
 235             {
 236                 delete[] pdblInI[i][j];
 237             }
 238         }
 239
 240         if (pdblInI[i])
 241         {
 242             delete[] pdblInI[i];
 243         }
 244     }
 245
 246     /*** retrun output arguments ***/
 247     if (iError)
 248     {
 249         return types::Function::Error;
 250     }
 251
 252     out.push_back(pDblOut);
 253     return types::Function::OK;
 254 }
 255 /*--------------------------------------------------------------------------*/