scilab/modules/arnoldi/src/c/eigs.c

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2012 -Scilab Enterprises - Adeline CARNIS
   4 *
   5 * This file must be used under the terms of the CeCILL.
   6 * This source file is licensed as described in the file COPYING, which
   7 * you should have received as part of this distribution.  The terms
   8 * are also available at
   9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  10 *
  11 */
  12 /*--------------------------------------------------------------------------*/
  13 #include <string.h>
  14 #include <math.h>
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 #include "eigs.h"
  18 #include "stack-c.h"
  19 #include "MALLOC.h"
  20 #include "sciprint.h"
  21 #include "eigs_dependencies.h"
  22 /*--------------------------------------------------------------------------*/
  23
  24 /*--------------------------------------------------------------------------*/
  25 // dgemm performs one of the matrix-matrix operations
  26 extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k,
  27                       double* alpha, double* A, int* lda, double* B, int* ldb,
  28                       double* beta, double* C, int* ldc);
  29 // zgemm performs one of the matrix-matrix operations
  30 extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k,
  31                       doublecomplex* alpha, doublecomplex* A, int* lda, doublecomplex* B,
  32                       int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc);
  33
  34 // dgemv performs matrix-vector operations
  35 extern int C2F(dgemv) (char* trans, int* m, int* n, double* alpha, double* A, int* lda,
  36                        double* x, int* incx, double* beta, double* y, int* incy);
  37 extern int C2F(zgemv) (char* trans, int* m, int* n, doublecomplex* alpha, doublecomplex* A,
  38                        int* lda, doublecomplex* x, int* incx, doublecomplex* beta, doublecomplex* y, int* incy);
  39
  40 // dgetrf computes an LU factorization of a general M by N matrix A (double) using partial pivoting with row interchanges
  41 extern int C2F(dgetrf)(int* m, int* n, double* A, int* lda, int* ipiv, int* info);
  42
  43 // zgetrf computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges
  44 extern int C2F(zgetrf)(int* m, int* n, doublecomplex* A, int* lda, int* ipiv, int* info);
  45 // dgetrs solves a linear system using the factors computed by dgetrf
  46 extern int C2F(dgetrs) (char* trans, int* n, int* nrhs, double* A, int *lda, int* ipiv, double* B, int* ldb, int* info);
  47 // zgetrs solves a linear system using the factors computed by zgetrf
  48 extern int C2F(zgetrs) (char* trans, int* n, int* nrhs, doublecomplex* AC, int* lda, int* ipiv, doublecomplex* B, int* ldb, int* info);
  49
  50 // dpotrf computes the cholesky factorization of a real symmetric positive definite matrix A
  51 extern int C2F(dpotrf)(char* uplo, int* n, double* A, int* lda, int* info);
  52
  53 // zpotrf computes the cholesky factorization of a real hermitian positive definite matrix A
  54 extern int C2F(zpotrf)(char* uplo, int* n, doublecomplex* A, int* lda, int* info);
  55
  56 // dtrsm solves a triangular linear system
  57 extern int C2F(dtrsm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, double* alpha, double* A, int* lda, double* B, int* ldb);
  58 // ztrsm solve a triangular linear system
  59 extern int C2F(ztrsm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, doublecomplex* alpha, doublecomplex* A, int* lda, doublecomplex* B, int* ldb);
  60 // dsyrk does a rank k symmetric update
  61 extern int C2F(dsyrk) (char* uplo, char* trans, int* n, int* k, double* alpha,
  62                        double* A, int* lda, double* beta, double* B, int* ldb);
  63 // ztrmm multiply by a triangular matrix
  64 extern int C2F(ztrmm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, doublecomplex* alphac,
  65                        doublecomplex* A, int* lda, doublecomplex* B, int* ldb);
  66 // ztrmv multiply a vector by a triangular matrix
  67 extern int C2F(ztrmv) (char* uplo, char* trans, char* diag, int* n, doublecomplex* A, int* lda, doublecomplex* x, int* incx);
  68 // dtrmv multiply a vector by a triangular matrix
  69 extern int C2F(dtrmv) (char* uplo, char* trans, char* diag, int* n, double* A, int* lda, double* x, int* incx);
  70 /*--------------------------------------------------------------------------*/
  71
  72 /*--------------------------------------------------------------------------*/
  73 extern int C2F(dsaupd)(int *ido, char *bmat, int *n, char *which, int *nev,
  74                        double *tol, double *resid, int *ncv, double *v,
  75                        int *ldv, int *iparam, int *ipntr, double *workd,
  76                        double *workl, int *lworkl, int *info);
  77 /*--------------------------------------------------------------------------*/
  78
  79 /*--------------------------------------------------------------------------*/
  80 extern int C2F(dseupd)(int *rvec, char *howmny, int *select, double *d,
  81                        double *z, int *ldz, double *sigma, char *bmat,
  82                        int *n, char *which, int *nev , double *tol,
  83                        double *resid, int *ncv, double *v , int *ldv,
  84                        int *iparam, int *ipntr, double *workd, double *workl,
  85                        int *lworkl, int *info);
  86 /*--------------------------------------------------------------------------*/
  87
  88 /*--------------------------------------------------------------------------*/
  89 extern int C2F(dnaupd)(int *ido, char *bmat, int *n, char *which, int *nev,
  90                        double *tol, double *resid, int *ncv, double *v,
  91                        int *ldv, int *iparam, int *ipntr, double *workd,
  92                        double *workl, int *lworkl, int *info);
  93 /*--------------------------------------------------------------------------*/
  94
  95 /*--------------------------------------------------------------------------*/
  96 extern int C2F(dneupd)(int *rvec, char *howmny, int *select, double *dr,
  97                        double *di, double *z, int *ldz, double *sigmar,
  98                        double *sigmai, double *workev, char *bmat, int *n,
  99                        char *which, int *nev, double *tol, double *resid,
 100                        int *ncv, double *v, int *ldv, int *iparam, int *ipntr,
 101                        double *workd, double *workl, int *lworkl, int *info);
 102 /*--------------------------------------------------------------------------*/
 103
 104 /*--------------------------------------------------------------------------*/
 105 extern int C2F(znaupd)(int * ido, char * bmat, int * n, char * which,
 106                        int * nev, double * tol, doublecomplex * resid,
 107                        int * ncv, doublecomplex * v, int * ldv, int * iparam,
 108                        int * ipntr, doublecomplex * workd,
 109                        doublecomplex * workl, int * lworkl, double * rwork,
 110                        int * info);
 111 /*--------------------------------------------------------------------------*/
 112
 113 /*--------------------------------------------------------------------------*/
 114 extern int C2F(zneupd)(int * rvec, char * howmny, int * select,
 115                        doublecomplex * d, doublecomplex * z, int * ldz,
 116                        doublecomplex * sigma, doublecomplex * workev,
 117                        char * bmat, int * n, char * which, int * nev,
 118                        double *  tol, doublecomplex * resid, int * ncv,
 119                        doublecomplex * v, int * ldv, int * iparam, int * ipntr,
 120                        doublecomplex * workd, doublecomplex * workl,
 121                        int * lworkl, double * rwork, int * info);
 122 /*--------------------------------------------------------------------------*/
 123
 124 static double alpha = 1.;
 125 static double beta = 0.;
 126
 127 static doublecomplex alphac = {1., 0.};
 128 static doublecomplex betac = {0., 0.};
 129
 130 int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym,
 131          double* B, doublecomplex* BC, int Bcomplex, int matB, int nev,
 132          doublecomplex SIGMA, char* which, double* maxiter, double* tol,
 133          double* NCV, double* RESID, doublecomplex* RESIDC, int* INFO,
 134          double* cholB, int INFO_EUPD, double* eigenvalue,
 135          double* eigenvector, doublecomplex* eigenvalueC,
 136          doublecomplex* eigenvectorC, int RVEC)
 137
 138 {
 139
 140     // GENERAL VARIABLES
 141     int i                       = 0;
 142     int j                       = 0;
 143     int k                       = 0;
 144     int l                       = 0;
 145     int INFO_CHOL       = 0;
 146     int INFO_LU         = 0;
 147     int INFO_INV    = 0;
 148     int iOne            = 1;
 149
 150     // VARIABLES DSAUPD, DNAUPD, ZNAUPD
 151     int LWORKL          = 0;
 152     int IDO                     = 0;
 153     int LDV                     = Max(1, N);
 154     int ncv                     = 0;
 155
 156     int IPARAM[11];
 157     int IPNTR[14];
 158
 159     double* V                   = NULL;
 160     doublecomplex* VC   = NULL;
 161
 162     double* WORKD                       = NULL;
 163     doublecomplex* WORKDC       = NULL;
 164
 165     double* WORKL                       = NULL;
 166     doublecomplex* WORKLC       = NULL;
 167
 168     double* RWORK                       = NULL;
 169
 170     char* bmat  = "I";
 171
 172     // VARIABLES DSEUPD, DNEUPD, ZNEUPD
 173     char* HOWMNY                = "A";
 174
 175     int* SELECT                 = NULL;
 176
 177     double* DI                  = NULL;
 178     double* DR                  = NULL;
 179     double* Z           = NULL;
 180
 181     double* WORKEV                      = NULL;
 182     doublecomplex* WORKEVC      = NULL;
 183
 184     doublecomplex mSIGMA = {.r = -SIGMA.r, .i = -SIGMA.i };
 185
 186     double* R         = NULL;
 187     doublecomplex* RC = NULL;
 188
 189     double* AMSB                        = NULL;
 190     doublecomplex* AMSBC        = NULL;
 191
 192     int* IPVT   = NULL;
 193
 194     double* temp = NULL;
 195     doublecomplex* tempC = NULL;
 196
 197     int oldnev = nev;
 198     int N2 = N * N;
 199
 200     IPARAM[0] = 1;
 201     IPARAM[2] = (int) maxiter[0];
 202     IPARAM[6] = 1; // by default mode = 1
 203
 204     // END VARIABLES
 205
 206     // MODE
 207     if (!strcmp(which, "SM") || (SIGMA.r != 0 || SIGMA.i != 0))
 208     {
 209         IPARAM[6] = 3;
 210         which = "LM";
 211     }
 212
 213     // BMAT
 214     if ((matB == 0) || (IPARAM[6] == 1))    // if B = [] or mode = 1 -> bmat = 'I' : standard eigenvalue problem
 215     {
 216         bmat = "I";
 217     }
 218     else   // generalized eigenvalue problem
 219     {
 220         bmat = "G";
 221     }
 222
 223     // NCV
 224     if (NCV == NULL)
 225     {
 226         if (Asym == 0 && !Acomplex && !Bcomplex) // if dnaupd  ncv = 2*nev+1
 227         {
 228             ncv = Max(2 * nev + 1, 20);
 229         }
 230         else // if dsaupd or znaupd ncv = 2*nev
 231         {
 232             ncv = Max(2 * nev, 20);
 233         }
 234         if (ncv > N)
 235         {
 236             ncv = N;
 237         }
 238     }
 239     else
 240     {
 241         ncv = (int) NCV[0];
 242         if (ncv <= nev || ncv > N) // Error
 243         {
 244             return -1;
 245         }
 246     }
 247
 248     // NEV
 249     if ((!Acomplex && !Bcomplex && Asym == 1 && nev >= N) || ((Acomplex || Bcomplex || !Asym) && nev >= N - 1))
 250     {
 251         return -2;
 252     }
 253
 254     if (matB != 0)
 255     {
 256         if (cholB[0]) // we already have the cholesky decomposition
 257         {
 258             R = B;
 259             RC = BC;
 260         }
 261         else
 262         {
 263             if (IPARAM[6] == 1)
 264             {
 265                 if (!Bcomplex) // B is real
 266                 {
 267                     R = (double *)malloc(N * N * sizeof(double));
 268                     memcpy(R, B, N * N * sizeof(double));
 269                     C2F(dpotrf) ("u", &N, R, &N, &INFO_CHOL);   // Compute the upper triangular matrix R
 270                     if (INFO_CHOL != 0) // Errors
 271                     {
 272                         free(R);
 273                         return -3;
 274                     }
 275                 }
 276                 else    // B is complex
 277                 {
 278                     RC = (doublecomplex *) malloc(N * N * sizeof(doublecomplex));
 279                     memcpy(RC, BC, N * N * sizeof(doublecomplex));
 280                     C2F(zpotrf) ("u", &N, RC, &N, &INFO_CHOL);  // Computes the upper triangular matrix
 281                     if (INFO_CHOL != 0)
 282                     {
 283                         free(RC);
 284                         return -3;
 285                     }
 286                 }
 287             }
 288         }
 289     }
 290
 291     // MAIN
 292     if (!Acomplex && !Bcomplex)         // A and B are not complex
 293     {
 294         if (IPARAM[6] == 3)     // if mode = 3
 295         {
 296             AMSB = (double*)malloc(N * N * sizeof(double));
 297             memcpy(AMSB, AR, N * N * sizeof(double));
 298             if (SIGMA.r != 0)
 299             {
 300                 // Compute LU decomposition AMSB = A - sigma*B
 301                 if (matB == 0) // if B = [] -> standard eigenvalue problem : A - sigma *I
 302                 {
 303                     for (i = 0 ; i < N ; i++)
 304                     {
 305                         AMSB[i + i * N] -= SIGMA.r;
 306                     }
 307                 }
 308                 else    // generalized eigenvalue problem
 309                 {
 310                     if (cholB[0])
 311                     {
 312                         C2F(dsyrk) ("u", "t", &N, &N, &mSIGMA.r, R, &N, &alpha, AMSB, &N);
 313                         if (!Asym)  //dsyrk does a symmetric update so we need to correct for the antisymmetric part
 314                         {
 315                             for (i = 0; i < N; i++)
 316                             {
 317                                 for (j = 0; j < i; j++)
 318                                 {
 319                                     AMSB[i + j * N] = AMSB[j + i * N] + AR[i + j * N] - AR[j + i * N];
 320                                 }
 321                             }
 322                         }
 323                     }
 324                     else
 325                     {
 326                         C2F(daxpy)(&N2, &mSIGMA.r, B, &iOne, AMSB, &iOne);
 327                     }
 328                 }
 329             }
 330
 331             // LU decomposition
 332             IPVT = (int*) calloc(N, sizeof(int));
 333             C2F(dgetrf)(&N, &N, AMSB, &N, IPVT, &INFO_LU);
 334             if (INFO_LU > 0)
 335             {
 336                 free(IPVT);
 337                 free(AMSB);
 338                 return -7;
 339             }
 340         }
 341
 342         if (Asym) // DSAUPD
 343         {
 344             LWORKL = ncv * ncv + 8 * ncv;
 345             WORKL = (double*) calloc(LWORKL, sizeof(double));
 346
 347         }
 348         else    // DNAUPD
 349         {
 350             LWORKL = 3 * ncv * (ncv + 2);
 351             WORKL = (double*) calloc(LWORKL, sizeof(double));
 352
 353         }
 354
 355         WORKD = (double*) calloc(3 * N, sizeof(double));
 356         V = (double*) calloc(N * ncv, sizeof(double));
 357
 358         if (IPARAM[6] == 1 && matB)
 359         {
 360             temp = (double*) malloc(N * sizeof(double));
 361         }
 362
 363         while (IDO != 99)
 364         {
 365             if (Asym) // DSAUPD
 366             {
 367                 C2F(dsaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0]);
 368             }
 369             else        // DNAUPD
 370             {
 371                 C2F(dnaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0]);
 372             }
 373
 374             if (INFO[0] == -1) //non critical error
 375             {
 376                 sciprint("%s: WARNING: Maximum number of iterations reached. Only %d eigenvalues converged.\n", "eigs", IPARAM[4]);
 377                 break;
 378             }
 379             else
 380             {
 381                 if (INFO[0] < 0)
 382                 {
 383                     if (R != B)
 384                     {
 385                         free(R);
 386                     }
 387                     free(IPVT);
 388                     free(AMSB);
 389                     free(WORKD);
 390                     free(WORKL);
 391                     free(V);
 392                     free(temp);
 393
 394                     return -4;
 395                 }
 396             }
 397
 398             if (IDO == -1 || IDO == 1 || IDO == 2)
 399             {
 400                 if (IPARAM[6] == 1) // mode = 1
 401                 {
 402                     if (IDO == 2)
 403                     {
 404                         memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
 405                     }
 406                     else        //IDO=1 or IDO=-1
 407                     {
 408                         if (matB == 0) // B = [] -> standard eigenvalue problem
 409                         {
 410                             // OP = A*x
 411                             if (Asym)
 412                             {
 413                                 C2F(dsymv) ("u", &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne);
 414                             }
 415                             else
 416                             {
 417                                 C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne);
 418                             }
 419                         }
 420                         else // generalized eigenvalue problem
 421                         {
 422                             // OP = inv(Rprime)*A*inv(R)*x
 423                             memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
 424                             C2F(dtrsm) ("l", "u", "n", "n", &N, &iOne, &alpha, R, &N, WORKD + IPNTR[1] - 1, &N);
 425                             memcpy(temp, WORKD + IPNTR[1] - 1, N * sizeof(double));
 426                             if (Asym)
 427                             {
 428                                 C2F(dsymv) ("u", &N, &alpha, AR, &N, temp, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne);
 429                             }
 430                             else
 431                             {
 432                                 C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, temp, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne);
 433                             }
 434                             C2F(dtrsm) ("l", "u", "t", "n", &N, &iOne, &alpha, R, &N, WORKD + IPNTR[1] - 1, &N);
 435                         }
 436                     }
 437                 }
 438                 else
 439                 {
 440                     if (IPARAM[6] == 3) // mode = 3
 441                     {
 442                         if (matB == 0) // B = [] -> standard eigenvalue problem
 443                         {
 444                             if (IDO == 2)
 445                             {
 446                                 // y = B*x where B = I so workd[ipntr[1]-1:ipntr[1]+N-1] = workd[ipntr[0]-1:ipntr[0]+N-1]
 447                                 memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
 448                             }
 449                             else
 450                             {
 451                                 // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1]
 452                                 memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
 453                                 C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV);
 454                             }
 455                         }
 456                         else  // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem
 457                         {
 458                             if (IDO == 2 || IDO == -1)
 459                             {
 460                                 if (cholB[0])   // workd[ipntr[1]-1:ipntr[1]+N-1] = Rprime * R * workd[ipntr[0]-1:ipntr[0]+N-1]
 461                                 {
 462                                     memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
 463                                     C2F(dtrmv) ("u", "n", "n", &N, B, &N, WORKD + IPNTR[1] - 1, &iOne);
 464                                     C2F(dtrmv) ("u", "t", "n", &N, B, &N, WORKD + IPNTR[1] - 1, &iOne);
 465                                 }
 466                                 else    //  workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1]
 467                                 {
 468                                     C2F(dgemv) ("n", &N, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &iOne, &beta, WORKD + IPNTR[1] - 1, &iOne);
 469                                 }
 470                             }
 471
 472                             if (IDO == -1)
 473                             {
 474                                 // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1]
 475                                 C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV);
 476                             }
 477                             else
 478                             {
 479                                 if (IDO == 1)
 480                                 {
 481                                     // computes workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1]
 482                                     memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[2] - 1, N * sizeof(double));
 483                                     C2F(dgetrs) ("n", &N, &iOne, AMSB, &N, IPVT, WORKD + IPNTR[1] - 1, &N, &INFO_INV);
 484                                 }
 485                             }
 486                         }
 487                     }
 488                     else
 489                     {
 490                         if (R != B)
 491                         {
 492                             free(R);
 493                         }
 494                         free(AMSB);
 495                         free(IPVT);
 496                         free(WORKD);
 497                         free(WORKL);
 498                         free(V);
 499                         free(temp);
 500
 501                         return -5;
 502                     }
 503                 }
 504             }
 505         } // END WHILE
 506         free(AMSB);
 507         free(IPVT);
 508         free(temp);
 509         SELECT = (int *)calloc(ncv, sizeof(int));
 510
 511         if (Asym) // DSEUPD
 512         {
 513             C2F(dseupd) (&RVEC, HOWMNY, SELECT, eigenvalue, eigenvector, &LDV,
 514                          &SIGMA.r, bmat, &N, which, &nev, tol, RESID, &ncv, V,
 515                          &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD);
 516
 517             if (INFO_EUPD != 0)
 518             {
 519                 if (R != B)
 520                 {
 521                     free(R);
 522                 }
 523                 free(WORKD);
 524                 free(WORKL);
 525                 free(V);
 526                 free(SELECT);
 527                 return -6;
 528             }
 529             else
 530             {
 531                 if (RVEC)
 532                 {
 533                     if (matB && IPARAM[6] == 1)
 534                     {
 535                         // we need to revert back to the original problem
 536                         // since we really solved for (y,\lambda) in R^{-T}Ay=\lambda y
 537                         //with y = Rx, so that x = R^{-1}y
 538                         C2F(dtrsm) ("l", "u", "n", "n", &N, &nev, &alpha, R, &N, eigenvector, &N);
 539                     }
 540                 }
 541             }
 542         }
 543         else    // DNEUPD
 544         {
 545             DR = (double *)calloc((nev + 1), sizeof(double));
 546             DI = (double *)calloc((nev + 1), sizeof(double));
 547             WORKEV = (double *)calloc(3 * ncv, sizeof(double));
 548
 549             RVEC = RVEC || (IPARAM[6] == 3 && SIGMA.i != 0);
 550
 551             if (RVEC)
 552             {
 553                 Z = (double *)calloc(N * (nev + 1), sizeof(double));
 554             }
 555
 556             C2F(dneupd) (&RVEC, HOWMNY, SELECT, DR, DI, Z, &LDV, &SIGMA.r,
 557                          &SIGMA.i, WORKEV, bmat, &N, which, &nev, tol, RESID,
 558                          &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD);
 559
 560             if (INFO_EUPD != 0)
 561             {
 562                 if (R != B)
 563                 {
 564                     free(R);
 565                 }
 566                 free(WORKD);
 567                 free(WORKL);
 568                 free(V);
 569                 free(DR);
 570                 free(DI);
 571                 free(Z);
 572                 free(WORKEV);
 573                 free(SELECT);
 574                 return -6;
 575             }
 576             else
 577             {
 578                 if (Z && matB && IPARAM[6] == 1)
 579                 {
 580                     // we need to revert back to the original problem
 581                     // since we really solved for (y,\lambda) in R^{-T}Ay=\lambda y
 582                     //with y = Rx, so that x = R^{-1}y
 583                     C2F(dtrsm) ("l", "u", "n", "n", &N, &nev, &alpha, R, &N, Z, &N);
 584                 }
 585                 //we use oldnev, because dneupd increases nev by one sometimes.
 586                 process_dneupd_data(DR, DI, Z, N, oldnev, AR, eigenvalueC,
 587                                     eigenvectorC, (IPARAM[6] == 3) && (SIGMA.i != 0));
 588
 589                 free(DR);
 590                 free(DI);
 591                 free(Z);
 592                 free(WORKEV);
 593             }
 594         }
 595
 596         free(V);
 597         free(WORKD);
 598         free(WORKL);
 599         free(SELECT);
 600         if (R != B)
 601         {
 602             free(R);
 603         }
 604     }
 605     else // A or/and B complex
 606     {
 607         if (IPARAM[6] == 3)     // mode = 3
 608         {
 609             AMSBC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
 610
 611             if (SIGMA.r != 0 || SIGMA.i != 0)
 612             {
 613                 if (matB == 0)  // standard eigenvalue problem
 614                 {
 615                     memcpy(AMSBC, AC, N * N * sizeof(doublecomplex));
 616                     for (i = 0 ; i < N ; i++)
 617                     {
 618                         AMSBC[i + i * N].r -= SIGMA.r;
 619                         AMSBC[i + i * N].i -= SIGMA.i;
 620                     }
 621                 }
 622                 else    // generalized eigenvalue problem
 623                 {
 624                     if (cholB[0])
 625                     {
 626                         memcpy(AMSBC, BC, N * N * sizeof(doublecomplex));
 627                         C2F(ztrmm)("l", "u", "c", "n", &N, &N, &mSIGMA, BC, &N, AMSBC, &N);
 628                         C2F(zaxpy)(&N2, &alphac, AC, &iOne, AMSBC, &iOne);
 629                     }
 630                     else
 631                     {
 632                         memcpy(AMSBC, AC, N * N * sizeof(doublecomplex));
 633                         C2F(zaxpy) (&N2, &mSIGMA, BC, &iOne, AMSBC, &iOne);
 634                     }
 635                 }
 636             }
 637             else
 638             {
 639                 memcpy(AMSBC, AC, N * N * sizeof(doublecomplex));
 640             }
 641
 642             // LU decomposition
 643             IPVT = (int*) calloc(N, sizeof(int));
 644             C2F(zgetrf) (&N, &N, AMSBC, &N, IPVT, &INFO_LU);
 645             if (INFO_LU > 0)
 646             {
 647                 free(IPVT);
 648                 free(AMSBC);
 649                 return(-7);
 650             }
 651         }
 652         LWORKL = 3 * ncv * ncv + 5 * ncv;
 653
 654         VC = (doublecomplex*) calloc(N * ncv, sizeof(doublecomplex));
 655         WORKLC = (doublecomplex*) calloc(LWORKL, sizeof(doublecomplex));
 656         WORKDC = (doublecomplex*) calloc(3 * N, sizeof(doublecomplex));
 657         RWORK = (double*) calloc(ncv, sizeof(double));
 658         if (IPARAM[6] == 1 && matB)
 659         {
 660             tempC = (doublecomplex*) malloc(N * sizeof(doublecomplex));
 661         }
 662
 663         while (IDO != 99)
 664         {
 665             C2F(znaupd)(&IDO, bmat, &N, which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, WORKLC, &LWORKL, RWORK, &INFO[0]);
 666
 667             if (INFO[0] == -1) //non critical error
 668             {
 669                 sciprint("%s: WARNING: Maximum number of iterations reached. Only %d eigenvalues converged.\n", "eigs", IPARAM[4]);
 670                 break;
 671             }
 672             else
 673             {
 674                 if (INFO[0] < 0)
 675                 {
 676                     if (RC != BC)
 677                     {
 678                         free(RC);
 679                     }
 680                     free(WORKDC);
 681                     free(WORKLC);
 682                     free(VC);
 683                     free(RWORK);
 684                     return -4;
 685                 }
 686             }
 687
 688             if (IDO == -1 || IDO == 1 || IDO == 2)
 689             {
 690                 if (IPARAM[6] == 1) // mode = 1
 691                 {
 692                     if (IDO == 2)
 693                     {
 694                         memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
 695                     }
 696                     else
 697                     {
 698                         if (matB == 0) // B = I
 699                         {
 700                             // OP = A*x
 701                             C2F(zgemv) ("n", &N, &N, &alphac, AC, &N, WORKDC + IPNTR[0] - 1, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne);
 702                         }
 703                         else
 704                         {
 705                             // OP = inv(RC')*A*inv(RC)*x
 706                             memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
 707                             C2F(ztrsm) ("l", "u", "n", "n", &N, &iOne, &alphac, RC, &N, WORKDC + IPNTR[1] - 1, &N);
 708                             memcpy(tempC, WORKDC + IPNTR[1] - 1, N * sizeof(doublecomplex));
 709                             C2F(zgemv) ("n", &N, &N, &alphac, AC, &N, tempC, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne);
 710                             C2F(ztrsm) ("l", "u", "c", "n", &N, &iOne, &alphac, RC, &N, WORKDC + IPNTR[1] - 1, &N);
 711                         }
 712                     }
 713                 }
 714                 else
 715                 {
 716                     if (IPARAM[6] == 3) // if mode = 3
 717                     {
 718                         if (matB == 0)  // B = [] -> matB is empty -> standard eigenvalue problem
 719                         {
 720                             if (IDO == 2)
 721                             {
 722                                 // y = B*x where B = I so workd[ipntr[1]-1:ipntr[1]+N-1] = workd[ipntr[0]-1:ipntr[0]+N-1]
 723                                 memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
 724                             }
 725                             else
 726                             {
 727                                 // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1]
 728                                 memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
 729                                 C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV);
 730                             }
 731
 732                         }
 733                         else  // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem
 734                         {
 735                             if (IDO == 2 || IDO == -1)
 736                             {
 737                                 if (cholB[0])   // workd[ipntr[1]-1:ipntr[1]+N-1] = RCprime * RC * workd[ipntr[0]-1:ipntr[0]+N-1]
 738                                 {
 739                                     memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
 740                                     C2F(ztrmv) ("u", "n", "n", &N, BC, &N, WORKDC + IPNTR[1] - 1, &iOne);
 741                                     C2F(ztrmv) ("u", "c", "n", &N, BC, &N, WORKDC + IPNTR[1] - 1, &iOne);
 742                                 }
 743                                 else    // workd[ipntr[1]-1:ipntr[1]+N-1] = B *workd[ipntr[0]-1:ipntr[0]+N-1]
 744                                 {
 745                                     C2F(zgemv) ("n", &N, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &iOne, &betac, WORKDC + IPNTR[1] - 1, &iOne);
 746                                 }
 747                             }
 748                             if (IDO == -1)
 749                             {
 750                                 // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1]
 751                                 C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV);
 752                             }
 753                             else
 754                             {
 755                                 if (IDO == 1)
 756                                 {
 757                                     /* compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1] */
 758                                     memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[2] - 1, N * sizeof(doublecomplex));
 759                                     C2F(zgetrs) ("n", &N, &iOne, AMSBC, &N, IPVT, WORKDC + IPNTR[1] - 1, &N, &INFO_INV);
 760                                 }
 761                             }
 762                         }       //END mode3
 763                     }
 764                     else
 765                     {
 766                         if (RC != BC)
 767                         {
 768                             free(RC);
 769                         }
 770                         free(WORKDC);
 771                         free(WORKLC);
 772                         free(VC);
 773                         free(RWORK);
 774                         free(tempC);
 775
 776                         return -5;
 777                     }
 778                 }
 779             }
 780         } // END WHILE
 781         free(tempC);
 782         free(IPVT);
 783         free(AMSBC);
 784
 785         SELECT = (int *)calloc(ncv, sizeof(int));
 786         WORKEVC = (doublecomplex *) calloc(2 * ncv, sizeof(doublecomplex));
 787
 788         C2F(zneupd) (&RVEC, HOWMNY, SELECT, eigenvalueC, eigenvectorC, &LDV, &SIGMA, WORKEVC, bmat, &N,
 789                      which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC,
 790                      WORKLC, &LWORKL, RWORK, &INFO_EUPD);
 791
 792         if (INFO_EUPD != 0)
 793         {
 794             if (RC != BC)
 795             {
 796                 free(RC);
 797             }
 798             free(WORKDC);
 799             free(WORKLC);
 800             free(VC);
 801             free(SELECT);
 802             free(WORKEVC);
 803             free(RWORK);
 804
 805             return -6;
 806         }
 807         else
 808         {
 809             if (RVEC)
 810             {
 811                 if (matB && IPARAM[6] == 1)
 812                 {
 813                     C2F(ztrsm) ("l", "u", "n", "n", &N, &nev, &alphac, RC, &N, eigenvectorC, &N);
 814                 }
 815             }
 816         }
 817
 818         free(SELECT);
 819         free(WORKEVC);
 820
 821         free(VC);
 822         free(WORKDC);
 823         free(WORKLC);
 824         free(RWORK);
 825         if (RC != BC)
 826         {
 827             free(RC);
 828         }
 829     }
 830
 831     return 0;
 832 }