1e78840b00a17fc9e524bd0719ca6d44e780b80d
[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, double* alpha, double* A, int* lda, double* B, int* ldb, double* beta, double* C, int* ldc);
27
28 // zgemm performs one of the matrix-matrix operations
29 extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k, doublecomplex* alpha, doublecomplex* A, int* lda,
30                       doublecomplex* B, int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc);
31
32 // dgetrf computes an LU factorization of a general M by N matrix A (double) using partial pivoting with row interchanges
33 extern int C2F(dgetrf)(int* m, int* n, double* A, int* lda, int* ipiv, int* info);
34
35 // zgetrd computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges
36 extern int C2F(zgetrf)(int* m, int* n, doublecomplex* A, int* lda, int* ipiv, int* info);
37
38 // dlaswp performs a series of row interchanges on the matrix A
39 // one row interchange is initiated for each of rows k1 through k2 of A
40 extern int C2F(dlaswp)(int* n, double* A, int* lda, int* k1, int* k2, int* ipiv, int* incx);
41
42 // dpotrf computes the cholesky factorization of a real symmetric positive definite matrix A
43 extern int C2F(dpotrf)(char* uplo, int* n, double* A, int* lda, int* info);
44
45 // zpotrf computes the cholesky factorization of a real hermitian positive definite matrix A
46 extern int C2F(zpotrf)(char* uplo, int* n, doublecomplex* A, int* lda, int* info);
47
48 // dgetri computes the inverse of a matrix using the LU factorization computed by dgetrf
49 extern int C2F(dgetri)(int* n, double* A, int* lda, int* ipiv, double* work, int* lworkl, int* info);
50
51 // zgetri computes the inverse of a matrix using the LU factorization computed by zgetrf
52 extern int C2F(zgetri)(int* n, doublecomplex* A, int* lda, int* ipiv, doublecomplex* work, int* lworkl, int* info);
53 /*--------------------------------------------------------------------------*/
54
55 /*--------------------------------------------------------------------------*/
56 extern int C2F(dsaupd)(int *ido, char *bmat, int *n, char *which, int *nev,
57                        double *tol, double *resid, int *ncv, double *v,
58                        int *ldv, int *iparam, int *ipntr, double *workd,
59                        double *workl, int *lworkl, int *info);
60 /*--------------------------------------------------------------------------*/
61
62 /*--------------------------------------------------------------------------*/
63 extern int C2F(dseupd)(int *rvec, char *howmny, int *select, double *d,
64                        double *z, int *ldz, double *sigma, char *bmat,
65                        int *n, char *which, int *nev , double *tol,
66                        double *resid, int *ncv, double *v , int *ldv,
67                        int *iparam, int *ipntr, double *workd, double *workl,
68                        int *lworkl, int *info, unsigned long rvec_length,
69                        unsigned long howmany_length,
70                        unsigned long bmat_length, unsigned long which_len);
71 /*--------------------------------------------------------------------------*/
72
73 /*--------------------------------------------------------------------------*/
74 extern int C2F(dnaupd)(int *ido, char *bmat, int *n, char *which, int *nev,
75                        double *tol, double *resid, int *ncv, double *v,
76                        int *ldv, int *iparam, int *ipntr, double *workd,
77                        double *workl, int *lworkl, int *info,
78                        unsigned long bmat_len, unsigned long which_len);
79 /*--------------------------------------------------------------------------*/
80
81 /*--------------------------------------------------------------------------*/
82 extern int C2F(dneupd)(int *rvec, char *howmny, int *select, double *dr,
83                        double *di, double *z, int *ldz, double *sigmar,
84                        double *sigmai, double *workev, char *bmat, int *n,
85                        char *which, int *nev, double *tol, double *resid,
86                        int *ncv, double *v, int *ldv, int *iparam, int *ipntr,
87                        double *workd, double *workl, int *lworkl, int *info);
88 /*--------------------------------------------------------------------------*/
89
90 /*--------------------------------------------------------------------------*/
91 extern int C2F(znaupd)(int * ido, char * bmat, int * n, char * which,
92                        int * nev, double * tol, doublecomplex * resid,
93                        int * ncv, doublecomplex * v, int * ldv, int * iparam,
94                        int * ipntr, doublecomplex * workd,
95                        doublecomplex * workl, int * lworkl, double * rwork,
96                        int * info);
97 /*--------------------------------------------------------------------------*/
98
99 /*--------------------------------------------------------------------------*/
100 extern int C2F(zneupd)(int * rvec, char * howmny, int * select,
101                        doublecomplex * d, doublecomplex * z, int * ldz,
102                        doublecomplex * sigma, doublecomplex * workev,
103                        char * bmat, int * n, char * which, int * nev,
104                        double *  tol, doublecomplex * resid, int * ncv,
105                        doublecomplex * v, int * ldv, int * iparam, int * ipntr,
106                        doublecomplex * workd, doublecomplex * workl,
107                        int * lworkl, double * rwork, int * info);
108 /*--------------------------------------------------------------------------*/
109
110 static double alpha = 1.;
111 static double beta = 0.;
112
113 static doublecomplex alphac = {1., 0.};
114 static doublecomplex betac = {0., 0.};
115
116 int eigs(double *AR, doublecomplex *AC, int N, int Acomplex, int Asym, double* B,
117          doublecomplex* BC, int Bcomplex, int matB, int nev, doublecomplex* SIGMA,
118          char* which, double* maxiter, double* tol, double* NCV, double* RESID, doublecomplex* RESIDC,
119          int* INFO, double* cholB, int INFO_EUPD, doublecomplex* eigenvalue, doublecomplex* eigenvector)
120 {
121
122     int index = 0;
123     // GENERAL VARIABLES
124     int i                       = 0;
125     int j                       = 0;
126     int k                       = 0;
127     int l                       = 0;
128     int INFO_CHOL       = 0;
129     int INFO_LU         = 0;
130     int k1                      = 1;
131     int iOne            = 1;
132
133     // VARIABLES DSAUPD, DNAUPD, ZNAUPD
134     int LWORKL          = 0;
135     int IDO                     = 0;
136     int LDV                     = Max(1, N);
137     int ncv                     = 0;
138
139     int* IPARAM         = NULL;
140     int* IPNTR          = NULL;
141
142     double* V                   = NULL;
143     doublecomplex* VC   = NULL;
144
145     double* WORKD                       = NULL;
146     doublecomplex* WORKDC       = NULL;
147
148     double* WORKL                       = NULL;
149     doublecomplex* WORKLC       = NULL;
150
151     double* RWORK                       = NULL;
152
153     char* bmat  = "I";
154
155     // VARIABLES DSEUPD, DNEUPD, ZNEUPD
156     int RVEC                    = 0;    // compute eigenvalues if RVEC = 1 also compute eigenvalues and eigenvectors
157     char* HOWMNY                = "A";
158
159     int* SELECT                 = NULL;
160
161     double* D                   = NULL;
162     double* DI                  = NULL;
163     double* DR                  = NULL;
164     doublecomplex* DC   = NULL;
165
166     double* WORKEV                      = NULL;
167     doublecomplex* WORKEVC      = NULL;
168
169     double* Z                   = NULL;
170     doublecomplex* ZC   = NULL;
171
172     double SIGMAR               = SIGMA[0].r;
173     double SIGMAI               = SIGMA[0].i;
174
175     double* R = (double*)malloc(N * N * sizeof(double));
176     double* Rprime = (double*)malloc(N * N * sizeof(double));
177
178     double* AMSB                        = NULL;
179     doublecomplex* AMSBC        = NULL;
180
181     double* L   = NULL;
182     double* U   = NULL;
183     double* E   = NULL;
184
185     doublecomplex* RC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
186     doublecomplex* RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
187
188     doublecomplex* LC   = NULL;
189     doublecomplex* UC   = NULL;
190     doublecomplex* EC   = NULL;
191
192     double* tmp_WORKD   = NULL;
193     doublecomplex* tmp_WORKDC   = NULL;
194
195     double* R_Rprime    = NULL;
196     double* invR_A_invRprime    = NULL;
197     double* invU_invL_E = NULL;
198
199     doublecomplex* RC_RCprime                   = NULL;
200     doublecomplex* invRC_AC_invRCprime  = NULL;
201     doublecomplex* invUC_invLC_EC               = NULL;
202
203     int* IPVT   = NULL;
204
205     IPARAM = (int*)malloc(11 * sizeof(int));
206     memset(IPARAM, 0, 11 * sizeof(int));
207     IPARAM[0] = 1;
208     IPARAM[2] = (int) maxiter[0];
209     IPARAM[6] = 1; // by default mode = 1
210
211     IPNTR = (int*)malloc(14 * sizeof(int));
212     memset(IPNTR, 0, 14 * sizeof(int));
213
214     tmp_WORKD = (double*)malloc(N * sizeof(double));
215     memset(tmp_WORKD, 0, N * sizeof(double));
216     tmp_WORKDC = (doublecomplex*)malloc(N * sizeof(doublecomplex));
217     memset(tmp_WORKDC, 0, N * sizeof(doublecomplex));
218
219     // END VARIABLES
220
221     // Info to compute eigenvalues and eigenvectors
222     if (eigenvector != NULL)
223     {
224         RVEC = 1;
225     }
226
227     // MODE
228     if (!strcmp(which, "SM") || (SIGMAR != 0 || SIGMAI != 0))
229     {
230         IPARAM[6] = 3;
231         which = "LM";
232     }
233
234     // BMAT
235     if ((matB == 0) || (IPARAM[6] == 1)) // if B = [] or mode = 1 -> bmat = 'I' : standart eigenvalue problem
236     {
237         bmat = "I";
238     }
239     else   // generalized eigenvalue problem
240     {
241         bmat = "G";
242     }
243
244     // NCV
245     if (NCV == NULL)
246     {
247         if (Asym == 0 && !Acomplex && !Bcomplex) // if dnaupd  ncv = 2*nev+1
248         {
249             ncv = Max(2 * nev + 1, 20);
250         }
251         else // if dsaupd or znaupd ncv = 2*nev
252         {
253             ncv = Max(2 * nev, 20);
254         }
255         if (ncv > N)
256         {
257             ncv = N;
258         }
259     }
260     else
261     {
262         ncv = (int) NCV[0];
263         if (ncv <= nev || ncv > N) // Error
264         {
265             free(IPARAM);
266             free(IPNTR);
267             free(R);
268             free(Rprime);
269             free(RC);
270             free(RCprime);
271             free(tmp_WORKD);
272             free(tmp_WORKDC);
273             return -1;
274
275         }
276     }
277
278     // NEV
279     if ((!Acomplex && !Bcomplex && Asym == 1 && nev >= N) || (!(!Acomplex && !Bcomplex && Asym == 1) && nev >= N - 1))
280     {
281         free(IPARAM);
282         free(IPNTR);
283         free(R);
284         free(Rprime);
285         free(RC);
286         free(RCprime);
287         free(tmp_WORKD);
288         free(tmp_WORKDC);
289         return -2;
290     }
291
292     // B must be symmetric (hermitian) positive (semi-) positive
293     if (matB != 0)
294     {
295         if (cholB[0]) // Comparison between B and upper triangular matrix
296         {
297             if (!Bcomplex) // B is real
298             {
299                 for (i = 0 ; i < N ; i++)
300                 {
301                     for (j = i + 1 ; j < N ; j++)
302                     {
303                         if (B[j + i * N] != 0)
304                         {
305                             free(IPARAM);
306                             free(IPNTR);
307                             free(R);
308                             free(Rprime);
309                             free(RC);
310                             free(RCprime);
311                             free(tmp_WORKD);
312                             free(tmp_WORKDC);
313                             free(IPVT);
314                             return -3;
315                         }
316                     }
317                 }
318
319                 memcpy(Rprime, B, N * N * sizeof(double));
320
321                 // Compute the lower triangular matrix
322                 memset(R, 0, N * N * sizeof(double));
323                 for (i = 0 ; i < N ; i++)
324                 {
325                     for (j = 0 ; j < N  ; j++)
326                     {
327                         R[i * N + j] = B[j * N + i];
328                     }
329                 }
330             }
331             else        // if B is complex
332             {
333                 for (i = 0 ; i < N ; i++)
334                 {
335                     for (j = i + 1 ; j < N ; j++)
336                     {
337                         if (BC[j + i * N].r != 0 || BC[j + i * N].i != 0)
338                         {
339                             free(IPARAM);
340                             free(IPNTR);
341                             free(R);
342                             free(Rprime);
343                             free(RC);
344                             free(RCprime);
345                             free(tmp_WORKD);
346                             free(tmp_WORKDC);
347                             return -3;
348                         }
349                     }
350                 }
351
352                 memcpy(RCprime, BC, N * N * sizeof(doublecomplex));
353
354                 // Computes the lower triangular matrix BC
355                 memset(RC, 0, N * N * sizeof(doublecomplex));
356                 for (i = 0 ; i < N ; i++)
357                 {
358                     for (j = 0 ; j < N ; j++)
359                     {
360                         RC[i * N + j].r = BC[j * N + i].r;
361                         RC[i * N + j].i = (-1) * BC[j * N + i].i;
362                     }
363                 }
364             }
365
366         }
367     }
368
369     if (!cholB[0] && IPARAM[6] == 1 && matB != 0)
370     {
371         if (!Bcomplex) // B is real
372         {
373             memcpy(R, B, N * N * sizeof(double));
374             memcpy(Rprime, B, N * N * sizeof(double));
375
376             C2F(dpotrf)("L", &N, R, &N, &INFO_CHOL); // Compute the lower triangular matrix R
377             if (INFO_CHOL != 0) // Errors
378             {
379                 free(IPARAM);
380                 free(IPNTR);
381                 free(R);
382                 free(Rprime);
383                 free(RC);
384                 free(RCprime);
385                 free(tmp_WORKD);
386                 free(tmp_WORKDC);
387                 return -3;
388             }
389
390             C2F(dpotrf)("U", &N, Rprime, &N, &INFO_CHOL);   // Compute the upper triangular matrix Rprime
391             if (INFO_CHOL != 0)
392             {
393                 free(IPARAM);
394                 free(IPNTR);
395                 free(R);
396                 free(Rprime);
397                 free(RC);
398                 free(RCprime);
399                 free(tmp_WORKD);
400                 free(tmp_WORKDC);
401                 return -3;
402             }
403
404             for (j = 0 ; j < N ; j++)
405             {
406                 for (i = 0 ; i < j ; i++)
407                 {
408                     R[i + j * N] = 0;
409                     Rprime[j + i * N] = 0;
410                 }
411             }
412         }
413         else    // B is complex
414         {
415             memcpy(RC, BC, N * N * sizeof(doublecomplex));
416             memcpy(RCprime, BC, N * N * sizeof(doublecomplex));
417
418             C2F(zpotrf)("L", &N, RC, &N, &INFO_CHOL); // Computes the lower triangular matrix
419             if (INFO_CHOL != 0)
420             {
421                 free(IPARAM);
422                 free(IPNTR);
423                 free(R);
424                 free(Rprime);
425                 free(RC);
426                 free(RCprime);
427                 free(tmp_WORKD);
428                 free(tmp_WORKDC);
429                 return -3;
430             }
431
432             C2F(zpotrf)("U", &N, RCprime, &N, &INFO_CHOL);      // Computes the upper triangular matrix
433             if (INFO_CHOL != 0)
434             {
435                 free(IPARAM);
436                 free(IPNTR);
437                 free(R);
438                 free(Rprime);
439                 free(RC);
440                 free(RCprime);
441                 free(tmp_WORKD);
442                 free(tmp_WORKDC);
443                 return -3;
444             }
445
446             for (j = 0 ; j < N ; j++)
447             {
448                 for (i = 0 ; i < j ; i++)
449                 {
450                     RC[i + j * N].r = 0;
451                     RC[i + j * N].i = 0;
452
453                     RCprime[j + i * N].r = 0;
454                     RCprime[j + i * N].i = 0;
455                 }
456             }
457         }
458     }
459
460     // MAIN
461     if (!Acomplex && !Bcomplex)         // A and B are not complex
462     {
463         if (IPARAM[6] == 3)     // if mode = 3
464         {
465             AMSB = (double*)malloc(N * N * sizeof(double));
466             memcpy(AMSB, AR, N * N * sizeof(double));
467
468             // Compute LU decomposition AMSB = A - sigma*B
469             if (matB == 0) // if B = [] -> standart eigenvalue problem : A - sigma
470             {
471                 for (i = 0 ; i < N ; i++)
472                 {
473                     AMSB[i + i * N] = AMSB[i + i * N] - SIGMAR;
474                 }
475             }
476             else        // generalized eigenvalue problem
477             {
478                 if (cholB[0])
479                 {
480                     if (R_Rprime == NULL)
481                     {
482                         R_Rprime = (double*)malloc(N * N * sizeof(double));
483                         RtimesRprime(R_Rprime, R, Rprime, N);
484                     }
485
486                     for (i = 0 ; i < N * N ; i++)
487                     {
488                         AMSB[i] = AR[i] - (SIGMAR * R_Rprime[i]);
489                     }
490                 }
491                 else
492                 {
493                     for (i = 0 ; i < N * N ; i++)
494                     {
495                         AMSB[i] = AR[i] - (SIGMAR * B[i]);
496                     }
497                 }
498             }
499
500             // LU decomposition
501             IPVT = (int*) malloc(N * sizeof(int));
502             memset(IPVT, 0, N * sizeof(int));
503             C2F(dgetrf)(&N, &N, AMSB, &N, IPVT, &INFO_LU);
504
505             // Computes the lower triangular matrix L
506             L = (double*)malloc(N * N * sizeof(double));
507             memset(L, 0, N * N * sizeof(double));
508
509             for (i = 0 ; i < N ; i++)
510             {
511                 for (j = 0 ; j < i ; j++)
512                 {
513                     L[i + j * N] = AMSB[i + j * N];
514                 }
515
516                 L[i + i * N] = 1;
517             }
518
519             // Computes the upper triangular matrix U
520             U = (double*)malloc(N * N * sizeof(double));
521             memset(U, 0, N * N * sizeof(double));
522
523             for (j = 0 ; j < N ; j++)
524             {
525                 for (i = 0 ; i <= j ; i++)
526                 {
527                     //if(i <= j)
528                     U[i + j * N] = AMSB[i + j * N];
529                 }
530             }
531
532             // Computes the permutation matrix E
533             E = (double*) malloc(N * N * sizeof(double));
534             memset(E, 0, N * N * sizeof(double));
535
536             for (i = 0 ; i < N ; i++)
537             {
538                 E[i * N + i] = 1;
539             }
540
541             C2F(dlaswp)(&N, E, &N, &k1, &N, IPVT, &k1);
542
543             free(AMSB);
544             free(IPVT);
545         }
546
547         if (Asym) // DSAUPD
548         {
549             LWORKL = ncv * ncv + 8 * ncv;
550
551             WORKL = (double*)malloc(LWORKL * sizeof(double));
552             memset(WORKL, 0, LWORKL * sizeof(double));
553         }
554         else    // DNAUPD
555         {
556             LWORKL = 3 * ncv * (ncv + 2);
557
558             WORKL = (double*)malloc(LWORKL * sizeof(double));
559             memset(WORKL, 0, LWORKL * sizeof(double));
560         }
561
562         WORKD = (double*)malloc(3 * N * sizeof(double));
563         memset(WORKD, 0, 3 * N * sizeof(double));
564
565         V = (double*)malloc(N * ncv * sizeof(double));
566         memset(V, 0, N * ncv * sizeof(double));
567
568         while (IDO != 99)
569         {
570             if (Asym) // DSAUPD
571             {
572                 C2F(dsaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0]);
573             }
574             else        // DNAUPD
575             {
576                 C2F(dnaupd)(&IDO, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO[0], 1L, 2L);
577             }
578
579             if (INFO[0] < 0)
580             {
581                 free(IPARAM);
582                 free(IPNTR);
583                 free(R);
584                 free(Rprime);
585                 free(RC);
586                 free(RCprime);
587                 free(tmp_WORKD);
588                 free(tmp_WORKDC);
589                 free(WORKD);
590                 free(WORKL);
591                 free(V);
592                 free(U);
593                 free(L);
594                 free(E);
595                 if (R_Rprime != NULL)
596                 {
597                     free(R_Rprime);
598                 }
599                 return -4;
600             }
601
602             if (IDO == -1 || IDO == 1 || IDO == 2)
603             {
604                 if (IPARAM[6] == 1) // mode = 1
605                 {
606                     if (matB == 0) // B = [] -> standart eigenvalue problem
607                     {
608                         // OP = A*x
609                         C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, AR, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
610                     }
611                     else // generalized eigenvalue problem
612                     {
613                         // OP = inv(Rprime)*A*inv(R)*x
614                         if (invR_A_invRprime == NULL)
615                         {
616                             invR_A_invRprime = (double*)malloc(N * N * sizeof(double));
617                             invR_times_A_times_invRprime(invR_A_invRprime, R, AR, Rprime, N);
618                         }
619
620                         C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invR_A_invRprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
621                     }
622                 }
623                 else
624                 {
625                     if (IPARAM[6] == 3) // mode = 3
626                     {
627                         if (matB == 0) // B = [] -> standart eigenvalue problem
628                         {
629                             if (IDO == 2)
630                             {
631                                 // y = B*x where B = I so workd[ipntr[1]-1:ipntr[1]+N-1] = workd[ipntr[0]-1:ipntr[0]+N-1]
632                                 memcpy(WORKD + IPNTR[1] - 1, WORKD + IPNTR[0] - 1, N * sizeof(double));
633                             }
634                             else
635                             {
636                                 // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1]
637
638                                 if (invU_invL_E == NULL)
639                                 {
640                                     invU_invL_E = (double*)malloc(N * N * sizeof(double));
641                                     invU_times_invL_times_E(invU_invL_E, U, L, E, N);
642                                 }
643
644                                 C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
645                             }
646                         }
647                         else  // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem
648                         {
649                             if (IDO == 2)
650                             {
651                                 if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = R * Rprime * workd[ipntr[0]-1:ipntr[0]+N-1]
652                                 {
653                                     if (R_Rprime == NULL)
654                                     {
655                                         R_Rprime = (double*)malloc(N * N * sizeof(double));
656                                         RtimesRprime(R_Rprime, R, Rprime, N);
657                                     }
658
659                                     C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, R_Rprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
660                                 }
661                                 else    //  workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1]
662                                 {
663                                     C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
664                                 }
665                             }
666                             else
667                             {
668                                 if (IDO == -1)
669                                 {
670                                     if (cholB[0])  // workd[ipntr[1]-1:ipntr[1]+N-1] = R * Rprime * workd[ipntr[0]-1:ipntr[0]+N-1]
671                                     {
672                                         if (R_Rprime == NULL)
673                                         {
674                                             R_Rprime = (double*)malloc(N * N * sizeof(double));
675                                             RtimesRprime(R_Rprime, R, Rprime, N);
676                                         }
677
678                                         C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, R_Rprime, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
679                                     }
680                                     else        // workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1]
681                                     {
682                                         C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, B, &N, WORKD + IPNTR[0] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
683                                     }
684                                     // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1]
685
686                                     if (invU_invL_E == NULL)
687                                     {
688                                         invU_invL_E = (double*)malloc(N * N * sizeof(double));
689                                         invU_times_invL_times_E(invU_invL_E, U, L, E, N);
690                                     }
691
692                                     C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[1] - 1, &N, &beta, tmp_WORKD, &N);
693                                     memcpy(WORKD + IPNTR[1] - 1, tmp_WORKD, N * sizeof(double));
694                                 }
695                                 else
696                                 {
697                                     if (IDO == 1)
698                                     {
699                                         // computes workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1]
700                                         if (invU_invL_E == NULL)
701                                         {
702                                             invU_invL_E = (double*)malloc(N * N * sizeof(double));
703                                             invU_times_invL_times_E(invU_invL_E, U, L, E, N);
704                                         }
705
706                                         C2F(dgemm)("n", "n", &N, &iOne, &N, &alpha, invU_invL_E, &N, WORKD + IPNTR[2] - 1, &N, &beta, WORKD + IPNTR[1] - 1, &N);
707                                     }
708                                 }
709                             }
710                         }
711                     }
712                     else
713                     {
714                         free(IPARAM);
715                         free(IPNTR);
716                         free(R);
717                         free(Rprime);
718                         free(RC);
719                         free(RCprime);
720                         free(tmp_WORKD);
721                         free(tmp_WORKDC);
722                         free(WORKD);
723                         free(WORKL);
724                         free(V);
725                         free(U);
726                         free(L);
727                         free(E);
728
729                         return -5;
730                     }
731                 }
732             }
733         } // END WHILE
734
735         free(L);
736         free(U);
737         free(E);
738
739         if (R_Rprime != NULL)
740         {
741             free(R_Rprime);
742         }
743
744         if (invR_A_invRprime != NULL)
745         {
746             free(invR_A_invRprime);
747         }
748
749         if (invU_invL_E != NULL)
750         {
751             free(invU_invL_E);
752         }
753
754         SELECT = (int*)malloc(ncv * sizeof(int));
755         memset(SELECT, 0, ncv * sizeof(int));
756
757         if (Asym) // DSEUPD
758         {
759             D = (double*)malloc(nev * sizeof(double));
760             memset(D, 0, nev * sizeof(double));
761
762             Z = (double*)malloc(N * nev * sizeof(double));
763             memset(Z, 0, N * nev * sizeof(double));
764
765             C2F(dseupd)(&RVEC, HOWMNY, SELECT, D, Z, &LDV, &SIGMAR, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD, 1L, 1L, 1L, 2L);
766
767             if (INFO_EUPD != 0)
768             {
769                 free(IPARAM);
770                 free(IPNTR);
771                 free(R);
772                 free(Rprime);
773                 free(RC);
774                 free(RCprime);
775                 free(tmp_WORKD);
776                 free(tmp_WORKDC);
777                 free(WORKD);
778                 free(WORKL);
779                 free(V);
780                 free(D);
781                 free(Z);
782                 free(SELECT);
783                 return -6;
784             }
785             else
786             {
787                 for (i = 0 ; i < nev ; i++)
788                 {
789                     eigenvalue[i].r = D[i];
790                 }
791
792                 if (RVEC)
793                 {
794                     for (i = 0 ; i < N * nev ; i++)
795                     {
796                         eigenvector[i].r = Z[i];
797                     }
798                 }
799             }
800
801             free(D);
802             free(Z);
803         }
804         else    // DNEUPD
805         {
806             DR = (double*)malloc((nev + 1) * sizeof(double));
807             memset(DR, 0, (nev + 1) * sizeof(double));
808
809             DI = (double*)malloc((nev + 1) * sizeof(double));
810             memset(DI, 0, (nev + 1) * sizeof(double));
811
812             Z = (double*) malloc(N * (nev + 1) * sizeof(double));
813             memset(Z, 0, N * (nev + 1) * sizeof(double));
814
815             WORKEV = (double*)malloc(3 * ncv * sizeof(double));
816             memset(WORKEV, 0, 3 * ncv * sizeof(double));
817
818             C2F(dneupd)(&RVEC, HOWMNY, SELECT, DR, DI, Z, &LDV, &SIGMAR, &SIGMAI, WORKEV, bmat, &N, which, &nev, tol, RESID, &ncv, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO_EUPD);
819
820             if (INFO_EUPD != 0)
821             {
822                 free(IPARAM);
823                 free(IPNTR);
824                 free(R);
825                 free(Rprime);
826                 free(RC);
827                 free(RCprime);
828                 free(tmp_WORKD);
829                 free(tmp_WORKDC);
830                 free(WORKD);
831                 free(WORKL);
832                 free(V);
833                 free(DR);
834                 free(DI);
835                 free(Z);
836                 free(WORKEV);
837                 free(SELECT);
838                 return -6;
839             }
840             else
841             {
842                 for (i = 0 ; i < nev ; i++)
843                 {
844                     eigenvalue[i].r = DR[i];
845                     eigenvalue[i].i = DI[i];
846                 }
847
848                 if (RVEC)
849                 {
850                     i = 0;
851                     while (i <= (nev - 2))
852                     {
853                         for (j = 0; j < N; j++)
854                         {
855                             eigenvector[i * N + j].r = Z[i * N + j];
856                             eigenvector[i * N + j].i = Z[(i + 1) * N + j];
857                             eigenvector[(i + 1)*N + j].r = Z[i * N + j];
858                             eigenvector[(i + 1)*N + j].i = -Z[(i + 1) * N + j];
859                         }
860                         i = i + 2;
861                     }
862                 }
863
864             }
865             free(DR);
866             free(DI);
867             free(Z);
868             free(WORKEV);
869         }
870
871         free(V);
872         free(WORKD);
873         free(WORKL);
874         free(SELECT);
875     }
876     else // A or/and B complex
877     {
878         if (IPARAM[6] == 3)     // mode = 3
879         {
880             AMSBC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
881             memcpy(AMSBC, AC, N * N * sizeof(doublecomplex));
882             if (matB == 0)      // standart eigenvalue problem
883             {
884                 for (i = 0 ; i < N ; i++)
885                 {
886                     AMSBC[i + i * N].r = AMSBC[i + i * N].r - SIGMAR;
887                     AMSBC[i + i * N].i = AMSBC[i + i * N].i - SIGMAI;
888                 }
889             }
890             else        // generalized eigenvalue problem
891             {
892                 if (cholB[0])
893                 {
894                     if (RC_RCprime == NULL)
895                     {
896                         RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
897                         RCtimesRCprime(RC_RCprime, RC, RCprime, N);
898                     }
899
900                     for (i = 0 ; i < N * N ; i++)
901                     {
902                         AMSBC[i].r = AMSBC[i].r - (SIGMAR * RC_RCprime[i].r + SIGMAI * RC_RCprime[i].i);
903                         AMSBC[i].i = AMSBC[i].i - (SIGMAR * RC_RCprime[i].i + SIGMAI * RC_RCprime[i].r);
904                     }
905                 }
906                 else
907                 {
908                     for (i = 0 ; i < N * N ; i++)
909                     {
910                         AMSBC[i].r = AMSBC[i].r - (SIGMA[0].r * BC[i].r);
911                         AMSBC[i].i = AMSBC[i].i - (SIGMA[0].i * BC[i].i);
912                     }
913                 }
914             }
915
916             // LU decomposition
917             IPVT = (int*) malloc(N * sizeof(int));
918             memset(IPVT, 0, N * sizeof(int));
919             C2F(zgetrf)(&N, &N, AMSBC, &N, IPVT, &INFO_LU);
920
921             // Computes the lower triangular matrix L
922             LC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
923             memset(LC, 0, N * N * sizeof(doublecomplex));
924             for (i = 0 ; i < N ; i++)
925             {
926                 for (j = 0 ; j < i ; j++)
927                 {
928                     LC[i + j * N].r = AMSBC[i + j * N].r;
929                     LC[i + j * N].i = AMSBC[i + j * N].i;
930                 }
931                 LC[i + i * N].r = 1;
932                 LC[i + i * N].i = 0;
933             }
934
935             // Computes the upper triangular matrix U
936
937             UC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
938             memset(UC, 0, N * N * sizeof(doublecomplex));
939             for (j = 0 ; j < N ; j++)
940             {
941                 for (i = 0 ; i <= j ; i++)
942                 {
943                     UC[i + j * N].r = AMSBC[i + j * N].r;
944                     UC[i + j * N].i = AMSBC[i + j * N].i;
945                 }
946             }
947
948             // Computes the permutation matrix E
949             EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
950             E = (double*)malloc(N * N * sizeof(double));
951             memset(E, 0, N * N * sizeof(double));
952
953             for (i = 0 ; i < N ; i++)
954             {
955                 E[i * N + i] = 1;
956             }
957
958             C2F(dlaswp)(&N, E, &N, &k1, &N, IPVT, &k1);
959
960             memset(EC, 0, N * N * sizeof(doublecomplex));
961             for (i = 0 ; i < N * N ; i++)
962             {
963                 EC[i].r = E[i];
964             }
965
966             free(AMSBC);
967             free(IPVT);
968         }
969
970         LWORKL = 3 * ncv * ncv + 5 * ncv;
971
972         VC = (doublecomplex*)malloc(N * ncv * sizeof(doublecomplex));
973         memset(VC, 0, N * ncv * sizeof(doublecomplex));
974
975         WORKLC = (doublecomplex*)malloc(LWORKL * sizeof(doublecomplex));
976         memset(WORKLC, 0, LWORKL * sizeof(doublecomplex));
977
978         WORKDC = (doublecomplex*)malloc(3 * N * sizeof(doublecomplex));
979         memset(WORKDC, 0, 3 * N * sizeof(doublecomplex));
980
981         RWORK = (double*)malloc(ncv * sizeof(double));
982         memset(RWORK, 0, ncv * sizeof(double));
983
984         while (IDO != 99)
985         {
986             C2F(znaupd)(&IDO, bmat, &N, which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, WORKLC, &LWORKL, RWORK, &INFO[0]);
987
988             if (INFO[0] < 0)
989             {
990                 free(IPARAM);
991                 free(IPNTR);
992                 free(R);
993                 free(Rprime);
994                 free(RC);
995                 free(RCprime);
996                 free(tmp_WORKD);
997                 free(tmp_WORKDC);
998                 free(LC);
999                 free(UC);
1000                 free(EC);
1001                 free(E);
1002                 free(WORKDC);
1003                 free(WORKLC);
1004                 free(VC);
1005                 free(RWORK);
1006                 if (RC_RCprime != NULL)
1007                 {
1008                     free(RC_RCprime);
1009                 }
1010                 return -4;
1011             }
1012
1013             if (IDO == -1 || IDO == 1 || IDO == 2)
1014             {
1015                 if (IPARAM[6] == 1) // mode = 1
1016                 {
1017                     if (matB == 0) // B = I
1018                     {
1019                         // OP = A*x
1020                         C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, AC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1021                     }
1022                     else
1023                     {
1024                         // OP = inv(R')*A*inv(R)*x
1025                         if (invRC_AC_invRCprime == NULL)
1026                         {
1027                             invRC_AC_invRCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1028                             invRC_times_AC_times_invRCprime(invRC_AC_invRCprime, RC, AC,  RCprime, N);
1029                         }
1030
1031                         C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invRC_AC_invRCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1032                     }
1033                 }
1034                 else
1035                 {
1036                     if (IPARAM[6] == 3) // si mode = 3
1037                     {
1038                         if (matB == 0)  // B = [] -> matB is empty -> standart eigenvalue problem
1039                         {
1040                             if (IDO == 2)
1041                             {
1042                                 // y = B*x where B = I so workd[ipntr[1]-1:ipntr[1]+N-1] = workd[ipntr[0]-1:ipntr[0]+N-1]
1043                                 memcpy(WORKDC + IPNTR[1] - 1, WORKDC + IPNTR[0] - 1, N * sizeof(doublecomplex));
1044                             }
1045                             else
1046                             {
1047                                 // workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[0]-1:ipntr[0]+N-1]
1048                                 if (invUC_invLC_EC == NULL)
1049                                 {
1050                                     invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1051                                     invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N);
1052                                 }
1053                                 C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1054                             }
1055
1056                         }
1057                         else  // matB == 1 so B is not empty and bmat = 'G'-> generalized eigenvalue problem
1058                         {
1059                             if (IDO == 2)
1060                             {
1061                                 if (cholB[0]) // workd[ipntr[1]-1:ipntr[1]+N-1] = RC * RCprime * workd[ipntr[0]-1:ipntr[0]+N-1]
1062                                 {
1063                                     if (RC_RCprime == NULL)
1064                                     {
1065                                         RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1066                                         RCtimesRCprime(RC_RCprime, RC, RCprime, N);
1067                                     }
1068
1069                                     C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, RC_RCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1070                                 }
1071                                 else    // workd[ipntr[1]-1:ipntr[1]+N-1] = B *workd[ipntr[0]-1:ipntr[0]+N-1]
1072                                 {
1073                                     C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1074                                 }
1075                             }
1076                             else
1077                             {
1078                                 if (IDO == -1)
1079                                 {
1080                                     if (cholB[0])  // workd[ipntr[1]-1:ipntr[1]+N-1] = RC*RCprime*workd[ipntr[0]-1:ipntr[0]+N-1]
1081                                     {
1082                                         if (RC_RCprime == NULL)
1083                                         {
1084                                             RC_RCprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1085                                             RCtimesRCprime(RC_RCprime, RC, RCprime, N);
1086                                         }
1087
1088                                         C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, RC_RCprime, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1089                                     }
1090                                     else        // workd[ipntr[1]-1:ipntr[1]+N-1] = B * workd[ipntr[0]-1:ipntr[0]+N-1]
1091                                     {
1092                                         C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, BC, &N, WORKDC + IPNTR[0] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1093                                     }
1094
1095                                     // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[1]-1:ipntr[1]+N-1]
1096                                     if (invUC_invLC_EC == NULL)
1097                                     {
1098                                         invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1099                                         invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N);
1100                                     }
1101
1102                                     C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[1] - 1, &N, &betac, tmp_WORKDC, &N);
1103                                     memcpy(WORKDC + IPNTR[1] - 1, tmp_WORKDC, N * sizeof(doublecomplex*));
1104                                 }
1105                                 else
1106                                 {
1107                                     if (IDO == 1)
1108                                     {
1109                                         // compute workd[ipntr[1]-1:ipntr[1]+N-1] = inv(U)*inv(L)*inv(P)*workd[ipntr[2]-1:ipntr[2]+N-1]
1110                                         if (invUC_invLC_EC == NULL)
1111                                         {
1112                                             invUC_invLC_EC = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
1113                                             invUC_times_invLC_times_EC(invUC_invLC_EC, UC, LC, EC, N);
1114                                         }
1115
1116                                         C2F(zgemm)("n", "n", &N, &iOne, &N, &alphac, invUC_invLC_EC, &N, WORKDC + IPNTR[2] - 1, &N, &betac, WORKDC + IPNTR[1] - 1, &N);
1117                                     }
1118                                 }
1119                             }
1120                         }
1121                     }
1122                     else
1123                     {
1124                         free(IPARAM);
1125                         free(IPNTR);
1126                         free(R);
1127                         free(Rprime);
1128                         free(RC);
1129                         free(RCprime);
1130                         free(tmp_WORKD);
1131                         free(tmp_WORKDC);
1132                         free(LC);
1133                         free(UC);
1134                         free(EC);
1135                         free(E);
1136                         free(WORKDC);
1137                         free(WORKLC);
1138                         free(VC);
1139                         free(RWORK);
1140                         if (RC_RCprime != NULL)
1141                         {
1142                             free(RC_RCprime);
1143                         }
1144
1145                         if (invRC_AC_invRCprime != NULL)
1146                         {
1147                             free(invRC_AC_invRCprime);
1148                         }
1149
1150                         if (invUC_invLC_EC != NULL)
1151                         {
1152                             free(invUC_invLC_EC);
1153                         }
1154                         return -5;
1155                     }
1156                 }
1157             }
1158         } // END WHILE
1159         free(LC);
1160         free(UC);
1161         free(EC);
1162         free(E);
1163
1164         if (RC_RCprime != NULL)
1165         {
1166             free(RC_RCprime);
1167         }
1168
1169         if (invRC_AC_invRCprime != NULL)
1170         {
1171             free(invRC_AC_invRCprime);
1172         }
1173
1174         if (invUC_invLC_EC != NULL)
1175         {
1176             free(invUC_invLC_EC);
1177         }
1178
1179         SELECT = (int*)malloc(ncv * sizeof(int));
1180         memset(SELECT, 0, ncv * sizeof(int));
1181
1182         DC = (doublecomplex*)malloc((nev + 1) * sizeof(doublecomplex));
1183         memset(DC, 0, (nev + 1) * sizeof(doublecomplex));
1184
1185         ZC = (doublecomplex*)malloc(N * nev * sizeof(doublecomplex));
1186         memset(ZC, 0, N * nev * sizeof(doublecomplex));
1187
1188         WORKEVC = (doublecomplex*)malloc(2 * ncv * sizeof(doublecomplex));
1189         memset(WORKEVC, 0, 2 * ncv * sizeof(doublecomplex));
1190
1191         C2F(zneupd)(&RVEC, HOWMNY, SELECT, DC, ZC, &LDV, SIGMA, WORKEVC, bmat, &N, which, &nev, tol, RESIDC, &ncv, VC, &LDV, IPARAM, IPNTR, WORKDC, WORKLC, &LWORKL, RWORK, &INFO_EUPD);
1192         if (INFO_EUPD != 0)
1193         {
1194             free(IPARAM);
1195             free(IPNTR);
1196             free(R);
1197             free(Rprime);
1198             free(RC);
1199             free(RCprime);
1200             free(tmp_WORKD);
1201             free(tmp_WORKDC);
1202             free(WORKDC);
1203             free(WORKLC);
1204             free(VC);
1205             free(SELECT);
1206             free(DC);
1207             free(ZC);
1208             free(WORKEVC);
1209             free(RWORK);
1210
1211             return -6;
1212         }
1213         else
1214         {
1215             if (!RVEC)
1216             {
1217                 for (i = 0; i < nev; i++)
1218                 {
1219                     eigenvalue[i].r = DC[i].r;
1220                     eigenvalue[i].i = DC[i].i;
1221                 }
1222             }
1223             else  // return eigenvalues and eigenvectors
1224             {
1225                 memcpy(eigenvalue, DC, nev * sizeof(doublecomplex));
1226                 memcpy(eigenvector, ZC, N * nev * sizeof(doublecomplex));
1227             }
1228         }
1229
1230         free(SELECT);
1231         free(DC);
1232         free(ZC);
1233         free(WORKEVC);
1234
1235         free(VC);
1236         free(tmp_WORKDC);
1237         free(tmp_WORKD);
1238         free(WORKDC);
1239         free(WORKLC);
1240         free(RWORK);
1241     }
1242
1243     free(IPARAM);
1244     free(IPNTR);
1245
1246     free(R);
1247     free(Rprime);
1248     free(RC);
1249     free(RCprime);
1250
1251     return 0;
1252 }
1253