scilab/modules/linear_algebra/sci_gateway/cpp/sci_spec.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 #include "overload.hxx"
  19 #include "execvisitor.hxx"
  20
  21 extern "C"
  22 {
  23 #include "localization.h"
  24 #include "Scierror.h"
  25 #include "sciprint.h"
  26 #include "eigen.h"
  27 #include "issymmetric.h"
  28 #include "vfinite.h"
  29 }
  30
  31 bool isNoZeroImag(types::Double* _pDbl);
  32 /*--------------------------------------------------------------------------*/
  33
  34 types::Function::ReturnValue sci_spec(types::typed_list &in, int _iRetCount, types::typed_list &out)
  35 {
  36     types::Double* pDblA    = NULL;
  37     types::Double* pDblB    = NULL;
  38     double* pDataA          = NULL;
  39     double* pDataB          = NULL;
  40     bool symmetric          = FALSE;
  41     int iRet                = 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"), "spec", 1, 2);
  46         return types::Function::Error;
  47     }
  48
  49     if (_iRetCount > 2 * in.size())
  50     {
  51         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "spec", 1, 2 * in.size());
  52         return types::Function::Error;
  53     }
  54
  55     if ((in[0]->isDouble() == false))
  56     {
  57         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_spec";
  58         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
  59     }
  60     pDblA = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
  61
  62     if (pDblA->getCols() != pDblA->getRows())
  63     {
  64         Scierror(20, _("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "spec", 1);
  65         return types::Function::Error;
  66     }
  67
  68     if ((pDblA->getRows() == -1) || (pDblA->getCols() == -1)) // manage eye case
  69     {
  70         Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "spec", 1);
  71         return types::Function::Error;
  72     }
  73
  74     if ((pDblA->getCols() == 0) || (pDblA->getRows() == 0)) // size null
  75     {
  76         out.push_back(types::Double::Empty());
  77         for (int i = 1; i < _iRetCount; i++)
  78         {
  79             out.push_back(types::Double::Empty());
  80         }
  81         return types::Function::OK;
  82     }
  83
  84     if (in.size() == 1)
  85     {
  86         types::Double* pDblEigenValues  = NULL;
  87         types::Double* pDblEigenVectors = NULL;
  88
  89         if (pDblA->isComplex())
  90         {
  91             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
  92             if (!pDataA)
  93             {
  94                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
  95                 return types::Function::Error;
  96             }
  97         }
  98         else
  99         {
 100             pDataA = pDblA->getReal();
 101         }
 102
 103         int totalSize = pDblA->getSize();
 104         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
 105         {
 106             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
 107             return types::Function::Error;
 108         }
 109
 110         symmetric = isSymmetric(pDblA->getReal(), pDblA->getImg(), pDblA->isComplex(), pDblA->getRows(), pDblA->getCols()) == 1;
 111         int eigenValuesCols = (_iRetCount == 1) ? 1 : pDblA->getCols();
 112
 113         if (symmetric)
 114         {
 115             pDblEigenValues = new types::Double(pDblA->getCols(), eigenValuesCols);
 116             if (_iRetCount == 2)
 117             {
 118                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), pDblA->isComplex());
 119             }
 120         }
 121         else
 122         {
 123             pDblEigenValues  = new types::Double(pDblA->getCols(), eigenValuesCols, true);
 124             if (_iRetCount == 2)
 125             {
 126                 pDblEigenVectors = new types::Double(pDblA->getCols(), pDblA->getCols(), true);
 127             }
 128         }
 129
 130         if (pDblA->isComplex())
 131         {
 132             if (symmetric)
 133             {
 134                 iRet = iEigen1ComplexSymmetricM((doublecomplex*)pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
 135
 136                 if (iRet < 0)
 137                 {
 138                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
 139                     return types::Function::Error;
 140                 }
 141
 142                 if (iRet > 0)
 143                 {
 144                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
 145                     return types::Function::Error;
 146                 }
 147
 148                 if (_iRetCount == 2)
 149                 {
 150                     vGetPointerFromDoubleComplex((doublecomplex*)pDataA, pDblA->getSize() , pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
 151                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
 152                     out.push_back(pDblEigenVectors);
 153                 }
 154                 out.push_back(pDblEigenValues);
 155             }
 156             else // not symmetric
 157             {
 158                 doublecomplex* pEigenValues = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
 159                 doublecomplex* pEigenVectors = pDblEigenVectors ? (doublecomplex*)MALLOC(sizeof(doublecomplex) * pDblA->getSize()) : NULL;
 160                 iRet = iEigen1ComplexM((doublecomplex*)pDataA, pDblA->getCols(), pEigenValues, pEigenVectors);
 161                 if (iRet < 0)
 162                 {
 163                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
 164                     return types::Function::Error;
 165                 }
 166
 167                 if (iRet > 0)
 168                 {
 169                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
 170                     return types::Function::Error;
 171                 }
 172
 173                 if (_iRetCount == 2)
 174                 {
 175                     expandZToDiagonalOfCMatrix(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
 176                     vGetPointerFromDoubleComplex(pEigenVectors, pDblA->getSize(), pDblEigenVectors->getReal(), pDblEigenVectors->getImg());
 177
 178                     FREE(pEigenVectors);
 179                     out.push_back(pDblEigenVectors);
 180                 }
 181                 else
 182                 {
 183                     vGetPointerFromDoubleComplex(pEigenValues, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg());
 184                 }
 185                 out.push_back(pDblEigenValues);
 186                 FREE(pEigenValues);
 187             }
 188         }
 189         else // real
 190         {
 191             if (symmetric)
 192             {
 193                 iRet = iEigen1RealSymmetricM(pDataA, pDblA->getCols(), (_iRetCount == 2), pDblEigenValues->getReal());
 194
 195                 if (iRet < 0)
 196                 {
 197                     Scierror(998, _("%s: On entry to ZGEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
 198                     return types::Function::Error;
 199                 }
 200
 201                 if (iRet > 0)
 202                 {
 203                     Scierror(24, _("%s: Convergence problem, %d off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n"), "spec", iRet);
 204                     return types::Function::Error;
 205                 }
 206
 207                 if (_iRetCount == 2)
 208                 {
 209                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
 210                     out.push_back(pDblA);
 211                 }
 212                 out.push_back(pDblEigenValues);
 213             }
 214             else // not symmetric
 215             {
 216                 iRet = iEigen1RealM(pDataA, pDblA->getCols(), pDblEigenValues->getReal(), pDblEigenValues->getImg(), pDblEigenVectors ? pDblEigenVectors->getReal() : NULL, pDblEigenVectors ? pDblEigenVectors->getImg() : NULL);
 217
 218                 if (iRet < 0)
 219                 {
 220                     Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
 221                     return types::Function::Error;
 222                 }
 223
 224                 if (iRet > 0)
 225                 {
 226                     Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of WR and WI contain eigenvalues which have converged.\n"), "spec", iRet);
 227                     return types::Function::Error;
 228                 }
 229
 230                 if (_iRetCount == 2)
 231                 {
 232                     expandToDiagonalOfMatrix(pDblEigenValues->getReal(), pDblA->getCols());
 233                     expandToDiagonalOfMatrix(pDblEigenValues->getImg(), pDblA->getCols());
 234                     out.push_back(pDblEigenVectors);
 235                 }
 236                 out.push_back(pDblEigenValues);
 237             }
 238         }
 239     }
 240
 241     if (in.size() == 2)
 242     {
 243         types::Double* pDblL            = NULL;
 244         types::Double* pDblR            = NULL;
 245         types::Double* pDblBeta         = NULL;
 246         types::Double* pDblAlpha        = NULL;
 247         doublecomplex* pL               = NULL;
 248         doublecomplex* pR               = NULL;
 249         doublecomplex* pBeta            = NULL;
 250         doublecomplex* pAlpha           = NULL;
 251         bool bIsComplex                 = false;
 252
 253         if ((in[1]->isDouble() == false))
 254         {
 255             std::wstring wstFuncName = L"%"  + in[1]->getShortTypeStr() + L"_spec";
 256             return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
 257         }
 258         pDblB = in[1]->getAs<types::Double>()->clone()->getAs<types::Double>();
 259         if ((pDblA->getRows() != pDblB->getRows()) && (pDblA->getCols() != pDblB->getCols()))
 260         {
 261             Scierror(999, _("%s: Arguments %d and %d must have equal dimensions.\n"), "spec", 1, 2);
 262             return types::Function::Error;
 263         }
 264
 265         //chekc if A and B are real complex or with imag part at 0
 266         if (isNoZeroImag(pDblA) == false && isNoZeroImag(pDblB) == false)
 267         {
 268             //view A and B as real matrix
 269             bIsComplex = false;
 270         }
 271         else
 272         {
 273             bIsComplex = pDblA->isComplex() || pDblB->isComplex();
 274         }
 275
 276         if (bIsComplex)
 277         {
 278             if (pDblA->isComplex() == false)
 279             {
 280                 pDblA->setComplex(true);
 281             }
 282
 283             if (pDblB->isComplex() == false)
 284             {
 285                 pDblB->setComplex(true);
 286             }
 287
 288             pDataA = (double*)oGetDoubleComplexFromPointer(pDblA->getReal(), pDblA->getImg(), pDblA->getSize());
 289             pDataB = (double*)oGetDoubleComplexFromPointer(pDblB->getReal(), pDblB->getImg(), pDblB->getSize());
 290
 291             if (!pDblA || !pDblB)
 292             {
 293                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "spec");
 294                 return types::Function::Error;
 295             }
 296         }
 297         else
 298         {
 299             pDataA = pDblA->getReal();
 300             pDataB = pDblB->getReal();
 301         }
 302
 303         int totalSize = pDblA->getSize();
 304
 305         if ((pDblA->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataA) : C2F(vfinite)(&totalSize, pDataA)) == false)
 306         {
 307             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 1);
 308             return types::Function::Error;
 309         }
 310
 311         if ((pDblB->isComplex() ? C2F(vfiniteComplex)(&totalSize, (doublecomplex*)pDataB) : C2F(vfinite)(&totalSize, pDataB)) == false)
 312         {
 313             Scierror(264, _("%s: Wrong value for input argument %d: Must not contain NaN or Inf.\n"), "spec", 2);
 314             return types::Function::Error;
 315         }
 316
 317         switch (_iRetCount)
 318         {
 319             case 4:
 320             {
 321                 pDblL = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
 322                 if (bIsComplex)
 323                 {
 324                     pL = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
 325                 }
 326             }
 327             case 3:
 328             {
 329                 pDblR = new types::Double(pDblA->getRows(), pDblA->getCols(), true);
 330                 if (bIsComplex)
 331                 {
 332                     pR = (doublecomplex*)MALLOC(pDblA->getSize() * sizeof(doublecomplex));
 333                 }
 334             }
 335             case 2:
 336             {
 337                 if (bIsComplex)
 338                 {
 339                     pBeta = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
 340                 }
 341                 pDblBeta = new types::Double(pDblA->getCols(), 1, pBeta ? true : false);
 342             }
 343             default : // case 1:
 344             {
 345                 if (bIsComplex)
 346                 {
 347                     pAlpha = (doublecomplex*)MALLOC(pDblA->getCols() * sizeof(doublecomplex));
 348                 }
 349                 pDblAlpha = new types::Double(pDblA->getCols(), 1, true);
 350             }
 351         }
 352
 353         if (bIsComplex)
 354         {
 355             iRet = iEigen2ComplexM((doublecomplex*)pDataA, (doublecomplex*)pDataB, pDblA->getCols(), pAlpha, pBeta, pR, pL);
 356         }
 357         else
 358         {
 359             iRet = iEigen2RealM(    pDataA, pDataB, pDblA->getCols(),
 360                                     pDblAlpha->getReal(), pDblAlpha->getImg(),
 361                                     pDblBeta ? pDblBeta->getReal()  : NULL,
 362                                     pDblR    ? pDblR->getReal()     : NULL,
 363                                     pDblR    ? pDblR->getImg()      : NULL,
 364                                     pDblL    ? pDblL->getReal()     : NULL,
 365                                     pDblL    ? pDblL->getImg()      : NULL);
 366         }
 367
 368         if (iRet > 0)
 369         {
 370             sciprint(_("Warning :\n"));
 371             sciprint(_("Non convergence in the QZ algorithm.\n"));
 372             sciprint(_("The top %d  x %d blocks may not be in generalized Schur form.\n"), iRet);
 373         }
 374
 375         if (iRet < 0)
 376         {
 377             Scierror(998, _("%s: On entry to ZHEEV parameter number  3 had an illegal value (lapack library problem).\n"), "spec", iRet);
 378             return types::Function::Error;
 379         }
 380
 381         if (iRet > 0)
 382         {
 383             if (bIsComplex)
 384             {
 385                 if (iRet <= pDblA->getCols())
 386                 {
 387                     Scierror(24, _("%s: The QZ iteration failed in DGGEV.\n"), "spec");
 388                 }
 389                 else
 390                 {
 391                     if (iRet ==  pDblA->getCols() + 1)
 392                     {
 393                         Scierror(999, _("%s: Other than QZ iteration failed in DHGEQZ.\n"), "spec");
 394                     }
 395                     if (iRet ==  pDblA->getCols() + 2)
 396                     {
 397                         Scierror(999, _("%s: Error return from DTGEVC.\n"), "spec");
 398                     }
 399                 }
 400             }
 401             else
 402             {
 403                 Scierror(24, _("%s: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements and %d+1:N of W contain eigenvalues which have converged.\n"), "spec", iRet);
 404             }
 405             return types::Function::Error;
 406         }
 407
 408         if (bIsComplex)
 409         {
 410             switch (_iRetCount)
 411             {
 412                 case 4:
 413                     vGetPointerFromDoubleComplex(pL, pDblA->getSize(), pDblL->getReal(), pDblL->getImg());
 414                 case 3:
 415                     vGetPointerFromDoubleComplex(pR, pDblA->getSize(), pDblR->getReal(), pDblR->getImg());
 416                 case 2:
 417                     vGetPointerFromDoubleComplex(pBeta, pDblA->getCols(), pDblBeta->getReal(), pDblBeta->getImg());
 418                 default : // case 1:
 419                     vGetPointerFromDoubleComplex(pAlpha, pDblA->getCols(), pDblAlpha->getReal(), pDblAlpha->getImg());
 420             }
 421         }
 422
 423         switch (_iRetCount)
 424         {
 425             case 1:
 426             {
 427                 out.push_back(pDblAlpha);
 428                 break;
 429             }
 430             case 2:
 431             {
 432                 out.push_back(pDblAlpha);
 433                 out.push_back(pDblBeta);
 434                 break;
 435             }
 436             case 3:
 437             {
 438                 out.push_back(pDblAlpha);
 439                 out.push_back(pDblBeta);
 440                 out.push_back(pDblR);
 441                 break;
 442             }
 443             case 4:
 444             {
 445                 out.push_back(pDblAlpha);
 446                 out.push_back(pDblBeta);
 447                 out.push_back(pDblL);
 448                 out.push_back(pDblR);
 449             }
 450         }
 451
 452         if (pAlpha)
 453         {
 454             vFreeDoubleComplexFromPointer(pAlpha);
 455         }
 456         if (pBeta)
 457         {
 458             vFreeDoubleComplexFromPointer(pBeta);
 459         }
 460         if (pL)
 461         {
 462             vFreeDoubleComplexFromPointer(pL);
 463         }
 464         if (pR)
 465         {
 466             vFreeDoubleComplexFromPointer(pR);
 467         }
 468         if (pDblB->isComplex())
 469         {
 470             vFreeDoubleComplexFromPointer((doublecomplex*)pDataB);
 471         }
 472
 473     } // if(in.size() == 2)
 474
 475     if (pDblA->isComplex())
 476     {
 477         vFreeDoubleComplexFromPointer((doublecomplex*)pDataA);
 478     }
 479
 480     return types::Function::OK;
 481 }
 482 /*--------------------------------------------------------------------------*/
 483 bool isNoZeroImag(types::Double* _pDbl)
 484 {
 485     double* pdbl = _pDbl->getImg();
 486     if (pdbl)
 487     {
 488         for (int i = 0 ; i < _pDbl->getSize() ; i++)
 489         {
 490             if (pdbl[i])
 491             {
 492                 return true;
 493             }
 494         }
 495     }
 496     return false;
 497 }
 498 /*--------------------------------------------------------------------------*/