scilab/modules/cacsd/sci_gateway/cpp/sci_ppol.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 extern "C"
  20 {
  21 #include "Scierror.h"
  22 #include "localization.h"
  23 #include "elem_common.h"
  24
  25     extern void C2F(ssxmc)( int*, int*, double*, int*, double*, int*, int*, int*, double*,
  26                             double*, double*, double*, double*, double*, int*);
  27
  28     extern void C2F(polmc)( int*, int*, int*, int*, double*, double*, double*, double*, double*, double*,
  29                             int*, int*, int*, double*, double*, double*, double*, double*, double*, double*);
  30
  31 }
  32
  33 /*--------------------------------------------------------------------------*/
  34 types::Function::ReturnValue sci_ppol(types::typed_list &in, int _iRetCount, types::typed_list &out)
  35 {
  36     types::Double* pDblA = NULL;
  37     types::Double* pDblB = NULL;
  38     types::Double* pDblP = NULL;
  39
  40     double dblEps = 0.1 * sqrt(nc_eps_machine());
  41
  42     int iOne   = 1;
  43     int iSizeP = 0;
  44     int iColB  = 0;
  45
  46     if (in.size() != 3)
  47     {
  48         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ppol", 3);
  49         return types::Function::Error;
  50     }
  51
  52     if (_iRetCount > 1)
  53     {
  54         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ppol", 1);
  55         return types::Function::Error;
  56     }
  57
  58     /*** get inputs arguments ***/
  59     // get poles
  60     if (in[2]->isDouble() == false)
  61     {
  62         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 3);
  63         return types::Function::Error;
  64     }
  65
  66     pDblP = in[2]->getAs<types::Double>();
  67     iSizeP = pDblP->getSize();
  68
  69     // get matrix B
  70     if (in[1]->isDouble() == false)
  71     {
  72         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 2);
  73         return types::Function::Error;
  74     }
  75
  76     pDblB = in[1]->getAs<types::Double>();
  77
  78     if (pDblB->isComplex())
  79     {
  80         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 2);
  81         return types::Function::Error;
  82     }
  83
  84     iColB = pDblB->getCols();
  85
  86     // get matrix A
  87     if (in[0]->isDouble() == false)
  88     {
  89         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ppol", 1);
  90         return types::Function::Error;
  91     }
  92
  93     pDblA = in[0]->getAs<types::Double>();
  94
  95     if (pDblA->isComplex())
  96     {
  97         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ppol", 1);
  98         return types::Function::Error;
  99     }
 100
 101     if (pDblA->getRows() != pDblA->getCols())
 102     {
 103         Scierror(999, _("%s: Wrong size for input argument #%d: A square matrix expected.\n"), "ppol", 1);
 104         return types::Function::Error;
 105     }
 106
 107     if (pDblA->getRows() != iSizeP || pDblB->getRows() != iSizeP)
 108     {
 109         Scierror(999, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "ppol");
 110         return types::Function::Error;
 111     }
 112
 113     if (iSizeP == 0)
 114     {
 115         out.push_back(types::Double::Empty());
 116         return types::Function::OK;
 117     }
 118
 119     /*** perform operations ***/
 120
 121     pDblA = pDblA->clone();
 122     pDblB = pDblB->clone();
 123     pDblP = pDblP->clone();
 124     pDblP->setComplex(true);
 125
 126     types::Double* pDblOut = new types::Double(iColB, iSizeP);
 127
 128     double* pdblG   = pDblOut->get();
 129     double* pdblZ   = new double[iSizeP * iSizeP];
 130     int iSizeWork   = std::max(iSizeP * iColB + 3 * iColB,
 131                                iColB * iColB + iColB * std::max(2, iColB) + 3 * iColB + 2 * iSizeP);
 132     double* pdblW   = new double[iSizeWork];
 133     int* piW        = new int[iSizeP];
 134
 135     double* pdblPReal = pDblP->get();
 136     double* pdblPImg  = pDblP->getImg();
 137
 138     int idc = 0;
 139     int inc = 0;
 140     // calcul de la forme canonique orthogonale
 141     C2F(ssxmc)(&iSizeP, &iColB, pDblA->get(), &iSizeP, pDblB->get(), &idc, &inc, piW, pdblZ,
 142                pdblW + iColB, pdblW + iColB + iSizeP * iColB, pdblW + 2 * iColB + iSizeP * iColB,
 143                pdblW, &dblEps, &iOne);
 144
 145     if (idc < iSizeP)
 146     {
 147         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
 148         delete[] pdblZ;
 149         delete[] pdblW;
 150         delete[] piW;
 151         pDblA->killMe();
 152         pDblB->killMe();
 153         pDblP->killMe();
 154         pDblOut->killMe();
 155         return types::Function::Error;
 156     }
 157
 158     double* pdblW1 = pdblW  + iColB + iColB * iColB;
 159     double* pdblW2 = pdblW1 + iColB * std::max(2, iColB);
 160     double* pdblW3 = pdblW2 + iSizeP;
 161     double* pdblW4 = pdblW3 + iSizeP;
 162     double* pdblW5 = pdblW4 + iColB;
 163
 164     // placement de poles
 165     int iErr = 0;
 166     C2F(polmc)( &iSizeP, &iColB, &iSizeP, &iColB, pDblA->get(), pDblB->get(), pdblG, pdblPReal, pdblPImg,
 167                 pdblZ, &inc, piW, &iErr, pdblW, pdblW + iColB, pdblW1, pdblW2, pdblW3, pdblW4, pdblW5);
 168
 169     // free memory
 170     delete[] pdblZ;
 171     delete[] pdblW;
 172     delete[] piW;
 173     pDblA->killMe();
 174     pDblB->killMe();
 175     pDblP->killMe();
 176
 177     if (iErr)
 178     {
 179         Scierror(999, _("%s: Uncontrollable system.\n"), "ppol");
 180         pDblOut->killMe();
 181         return types::Function::Error;
 182     }
 183
 184     /*** return output arguments ***/
 185     out.push_back(pDblOut);
 186     return types::Function::OK;
 187 }
 188 /*--------------------------------------------------------------------------*/