Bug #12138 fixed - eigs(A,B) returned incorrect eigenvectors for dense matrices.
[scilab.git] / 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 }