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                 int iSize1 = piRank[0][i] + 1;
 149                 int iSize2 = piRank[1][i] + 1;
 150                 int iSize3 = piRank[2][i] + 1;
 151                 C2F(residu)(pdblInR[0][i], &iSize1, pdblInR[1][i], &iSize2,
 152                             pdblInR[2][i], &iSize3, &v, &dblEps, &iErr);
 153                 if (iErr)
 154                 {
 155                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "residu");
 156                     throw iErr;
 157                 }
 158
 159                 pdblOut[i] = v;
 160             }
 161         }
 162         else
 163         {
 164             // complex case
 165             pDblOut = new types::Double(iRows[0], iCols[0], true);
 166             double* pdblOutR = pDblOut->get();
 167             double* pdblOutI = pDblOut->getImg();
 168
 169             for (int i = 0; i < 3; i++)
 170             {
 171                 if (pdblInI[i] == NULL)
 172                 {
 173                     pdblInI[i] = new double*[iSize];
 174                     for (int j = 0; j < iSize; j++)
 175                     {
 176                         int iLen = piRank[i][j] + 1;
 177                         pdblInI[i][j] = new double[iLen];
 178                         memset(pdblInI[i][j], 0x00, iLen * sizeof(double));
 179                     }
 180
 181                     isDeletable[i] = true;
 182                 }
 183             }
 184
 185             for (int i = 0; i < iRows[0] * iCols[0]; i++)
 186             {
 187                 int iErr    = 0;
 188                 double real = 0;
 189                 double imag = 0;
 190
 191                 C2F(wesidu)(pdblInR[0][i], pdblInI[0][i], (piRank[0]) + i,
 192                             pdblInR[1][i], pdblInI[1][i], (piRank[1]) + i,
 193                             pdblInR[2][i], pdblInI[2][i], (piRank[2]) + i,
 194                             &real, &imag, &dblEps, &iErr);
 195
 196                 if (iErr)
 197                 {
 198                     Scierror(78, _("%s: An error occured in '%s'.\n"), "residu", "wesidu");
 199                     throw iErr;
 200                 }
 201
 202                 pdblOutR[i] = real;
 203                 pdblOutI[i] = imag;
 204             }
 205         }
 206     }
 207     catch (int error)
 208     {
 209         iError = error;
 210     }
 211
 212     // free memory
 213     for (int i = 0; i < 3; i++)
 214     {
 215         if (pDblIn[i])
 216         {
 217             delete pDblIn[i];
 218         }
 219
 220         if (pPoly[i])
 221         {
 222             delete pPoly[i];
 223         }
 224
 225         if (piRank[i])
 226         {
 227             delete[] piRank[i];
 228         }
 229
 230         if (pdblInR[i])
 231         {
 232             delete[] pdblInR[i];
 233         }
 234
 235         if (isDeletable[i])
 236         {
 237             for (int j = 0; j < iSize; j++)
 238             {
 239                 delete[] pdblInI[i][j];
 240             }
 241         }
 242
 243         if (pdblInI[i])
 244         {
 245             delete[] pdblInI[i];
 246         }
 247     }
 248
 249     /*** retrun output arguments ***/
 250     if (iError)
 251     {
 252         return types::Function::Error;
 253     }
 254
 255     out.push_back(pDblOut);
 256     return types::Function::OK;
 257 }
 258 /*--------------------------------------------------------------------------*/