Bug #12138 fixed - eigs(A,B) returned incorrect eigenvectors for dense matrices.
[scilab.git] / scilab / modules / arnoldi / sci_gateway / c / sci_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 <math.h>
14 #include <string.h>
15 #include "stack-c.h"
16 #include "isanan.h"
17 #include "core_math.h"
18 #include "gw_arnoldi.h"
19 #include "localization.h"
20 #include "Scierror.h"
21 #include "api_scilab.h"
22 #include "stdio.h"
23 #include "stdlib.h"
24 #include "sciprint.h"
25 #include "doublecomplex.h"
26 #include "MALLOC.h"
27 #include "eigs.h"
28
29 int sci_eigs(char *fname, unsigned long fname_len)
30 {
31     int *piAddressVarOne        = NULL;
32     int iRowsOne                        = 0;
33     int iColsOne                        = 0;
34     double elemt1                       = 0;
35     double elemt2                       = 0;
36     double* Areal                       = NULL;
37     doublecomplex* Acplx        = NULL;
38     int Asym                            = 1;
39     int Acomplex                        = 0;
40     int N                                       = 0;
41
42     int *piAddressVarTwo        = NULL;
43     int iTypeVarTwo                     = 0;
44     int iRowsTwo                        = 0;
45     int iColsTwo                        = 0;
46     double* Breal                       = NULL;
47     doublecomplex* Bcplx        = NULL;
48     int Bcomplex                        = 0;
49     int matB                            = 0;
50
51     int *piAddressVarThree      = NULL;
52     double dblNEV                       = 0;
53     int iNEV                            = 0;
54
55     int *piAddressVarFour       = NULL;
56     int iTypeVarFour            = 0;
57     int iRowsFour                       = 0;
58     int iColsFour                       = 0;
59     int iLen                            = 0;
60     char* pstData                       = NULL;
61     doublecomplex SIGMA;
62
63     int *piAddressVarFive       = NULL;
64     double dblMAXITER           = 0;
65
66     int *piAddressVarSix        = NULL;
67     double dblTOL                       = 0;
68
69     int *piAddressVarSeven      = NULL;
70     int TypeVarSeven            = 0;
71     int RowsSeven                       = 0;
72     int ColsSeven                       = 0;
73     double* dblNCV                      = NULL;
74
75     int *piAddressVarEight      = NULL;
76     int iTypeVarEight       = 0;
77     double dblCHOLB                     = 0;
78     int iCHOLB              = 0;
79
80     int *piAddressVarNine       = NULL;
81     int iTypeVarNine            = 0;
82     int iRowsNine                       = 0;
83     int iColsNine                       = 0;
84     double* RESID                       = NULL;
85     doublecomplex* RESIDC       = NULL;
86
87     int *piAddressVarTen        = NULL;
88     int iINFO                           = 0;
89     int RVEC                = 0;
90     // Output arguments
91     double* eigenvalue      = NULL;
92     double* eigenvector     = NULL;
93     doublecomplex* eigenvalueC  = NULL;
94     doublecomplex* eigenvectorC = NULL;
95
96     double* mat_eigenvalue      = NULL;
97     doublecomplex* mat_eigenvalueC  = NULL;
98     int INFO_EUPD                                       = 0;
99     int error                                           = 0;
100
101     SciErr sciErr;
102     int iErr                            = 0;
103     int i                                       = 0;
104     int j                                       = 0;
105
106     CheckRhs(1, 10);
107     CheckLhs(0, 2);
108
109     /****************************************
110     *           First variable : A              *
111     *****************************************/
112
113     sciErr = getVarAddressFromPosition(pvApiCtx, 1, &piAddressVarOne);
114     if (sciErr.iErr)
115     {
116         printError(&sciErr, 0);
117         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 1);
118         return 0;
119     }
120
121     sciErr = getVarDimension(pvApiCtx, piAddressVarOne, &iRowsOne, &iColsOne);
122     //check if A is a square matrix
123     if (iRowsOne * iColsOne == 1 || iRowsOne != iColsOne)
124     {
125         Scierror(999, _("%s: Wrong type for input argument #%d: A square matrix expected.\n"), "eigs", 1);
126         return 0;
127     }
128
129     N = iRowsOne;
130
131     //check if A is complex
132     if (isVarComplex(pvApiCtx, piAddressVarOne))
133     {
134         sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddressVarOne, &iRowsOne, &iColsOne, &Acplx);
135         Acomplex = 1;
136     }
137     else
138     {
139         sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarOne, &iRowsOne, &iColsOne, &Areal);
140
141         for (i = 0; i < iColsOne; i++)
142         {
143             for (j = 0; j < i; j++)
144             {
145                 elemt1 = Areal[j + i * iColsOne];
146                 elemt2 = Areal[j * iColsOne + i];
147                 if (fabs(elemt1 - elemt2) > 0)
148                 {
149                     Asym = 0;
150                     break;
151                 }
152             }
153             if (Asym == 0)
154             {
155                 break;
156             }
157         }
158     }
159
160     /****************************************
161     *           Second variable : B             *
162     *****************************************/
163     sciErr = getVarAddressFromPosition(pvApiCtx, 2, &piAddressVarTwo);
164     if (sciErr.iErr)
165     {
166         printError(&sciErr, 0);
167         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 2);
168         return 0;
169     }
170
171     sciErr = getVarType(pvApiCtx, piAddressVarTwo, &iTypeVarTwo);
172     if (sciErr.iErr || iTypeVarTwo != sci_matrix)
173     {
174         printError(&sciErr, 0);
175         Scierror(999, _("%s: Wrong type for input argument #%d: An empty matrix or full or sparse square matrix expected.\n"), "eigs", 2);
176         return 0;
177     }
178
179     sciErr = getVarDimension(pvApiCtx, piAddressVarTwo, &iRowsTwo, &iColsTwo);
180     matB = iRowsTwo * iColsTwo;
181     if (matB && (iRowsTwo != iRowsOne || iColsTwo != iColsOne))
182     {
183         Scierror(999, _("%s: Wrong dimension for input argument #%d: B must have the same size as A.\n"), "eigs", 2);
184         return 0;
185     }
186
187     if (isVarComplex(pvApiCtx, piAddressVarTwo))
188     {
189         sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddressVarTwo, &iRowsTwo, &iColsTwo, &Bcplx);
190         Bcomplex = 1;
191     }
192     else
193     {
194         sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarTwo, &iRowsTwo, &iColsTwo, &Breal);
195     }
196
197     if (matB != 0)
198     {
199         if (Acomplex && !Bcomplex)
200         {
201             Bcplx = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
202             memset(Bcplx, 0, N * N * sizeof(doublecomplex));
203             Bcomplex = 1;
204             for (i = 0 ; i < N * N ;  i++)
205             {
206                 Bcplx[i].r = Breal[i];
207             }
208         }
209         if (!Acomplex && Bcomplex)
210         {
211             Acplx = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
212             memset(Acplx, 0, N * N * sizeof(doublecomplex));
213             Acomplex = 1;
214             for (i = 0 ; i < N * N ;  i++)
215             {
216                 Acplx[i].r = Areal[i];
217             }
218         }
219     }
220
221
222     /****************************************
223     *                           NEV                                     *
224     *****************************************/
225     sciErr = getVarAddressFromPosition(pvApiCtx, 3, &piAddressVarThree);
226     if (sciErr.iErr)
227     {
228         printError(&sciErr, 0);
229         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 3);
230         free(Acplx);
231         return 0;
232     }
233
234     iErr = getScalarDouble(pvApiCtx, piAddressVarThree, &dblNEV);
235     if (iErr)
236     {
237         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 3);
238         return 0;
239     }
240
241     if (isVarComplex(pvApiCtx, piAddressVarThree))
242     {
243         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 3);
244         return 0;
245     }
246
247     if (dblNEV != floor(dblNEV) || (dblNEV <= 0))
248     {
249         Scierror(999, _("%s: Wrong type for input argument #%d: k must be a positive integer.\n"), "eigs", 3);
250         return 0;
251     }
252
253     if (!finite(dblNEV))
254     {
255         Scierror(999, _("%s: Wrong value for input argument #%d: k must be in the range 1 to N.\n"), "eigs", 3);
256         return 0;
257     }
258
259
260     iNEV = (int)dblNEV;
261
262     /****************************************
263     *                   SIGMA AND WHICH                         *
264     *****************************************/
265     sciErr = getVarAddressFromPosition(pvApiCtx, 4, &piAddressVarFour);
266     if (sciErr.iErr)
267     {
268         printError(&sciErr, 0);
269         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 4);
270         return 0;
271     }
272
273     sciErr = getVarType(pvApiCtx, piAddressVarFour, &iTypeVarFour);
274     if (sciErr.iErr || (iTypeVarFour != sci_matrix && iTypeVarFour != sci_strings))
275     {
276         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 4);
277         return 0;
278     }
279
280     if (iTypeVarFour == sci_strings)
281     {
282         int iErr = getAllocatedSingleString(pvApiCtx, piAddressVarFour, &pstData);
283         if (iErr)
284         {
285             return 0;
286         }
287
288         if (strcmp(pstData, "LM") != 0 && strcmp(pstData, "SM") != 0  && strcmp(pstData, "LR") != 0 && strcmp(pstData, "SR") != 0 && strcmp(pstData, "LI") != 0
289                 && strcmp(pstData, "SI") != 0 && strcmp(pstData, "LA") != 0 && strcmp(pstData, "SA") != 0 && strcmp(pstData, "BE") != 0)
290         {
291             if (!Acomplex && Asym)
292             {
293                 Scierror(999, _("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of '%s', '%s', '%s', '%s' or '%s'.\n" ),
294                          "eigs", 4, "LM", "SM", "LA", "SA", "BE");
295                 return 0;
296             }
297             else
298             {
299                 Scierror(999, _("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of '%s', '%s', '%s', '%s', '%s' or '%s'.\n " ),
300                          "eigs", 4, "LM", "SM", "LR", "SR", "LI", "SI");
301                 return 0;
302             }
303         }
304
305         if ((Acomplex || !Asym) && (strcmp(pstData, "LA") == 0 || strcmp(pstData, "SA") == 0 || strcmp(pstData, "BE") == 0))
306         {
307             Scierror(999, _("%s: Invalid sigma value for complex or non symmetric problem.\n"), "eigs", 4);
308             return 0;
309         }
310
311         if (!Acomplex && Asym && (strcmp(pstData, "LR") == 0 || strcmp(pstData, "SR") == 0 || strcmp(pstData, "LI") == 0 || strcmp(pstData, "SI") == 0))
312         {
313             Scierror(999, _("%s: Invalid sigma value for real symmetric problem.\n"), "eigs", 4);
314             return 0;
315         }
316
317         SIGMA.r = 0;
318         SIGMA.i = 0;
319     }
320
321     if (iTypeVarFour == sci_matrix)
322     {
323         sciErr = getVarDimension(pvApiCtx, piAddressVarFour, &iRowsFour, &iColsFour);
324         if (iRowsFour * iColsFour != 1)
325         {
326             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 4);
327             return 0;
328         }
329
330         if (getScalarComplexDouble(pvApiCtx, piAddressVarFour, &SIGMA.r, &SIGMA.i))
331         {
332             printError(&sciErr, 0);
333             Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 4);
334             return 0;
335         }
336
337         if (C2F(isanan)(&SIGMA.r) || C2F(isanan)(&SIGMA.i))
338         {
339             Scierror(999, _("%s: Wrong type for input argument #%d: sigma must be a real.\n"), "eigs", 4);
340             return 0;
341         }
342
343         pstData = "LM";
344     }
345
346     /****************************************
347     *                           MAXITER                                 *
348     *****************************************/
349     sciErr = getVarAddressFromPosition(pvApiCtx, 5, &piAddressVarFive);
350     if (sciErr.iErr)
351     {
352         printError(&sciErr, 0);
353         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 5);
354         return 0;
355     }
356
357     iErr = getScalarDouble(pvApiCtx, piAddressVarFive, &dblMAXITER);
358     if (iErr)
359     {
360         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be a scalar.\n"), "eigs", 5, "opts.maxiter");
361         return 0;
362     }
363
364     if ((dblMAXITER != floor(dblMAXITER)) || (dblMAXITER <= 0))
365     {
366         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer positive value.\n"), "eigs", 5, "opts.maxiter");
367         return 0;
368     }
369
370     /****************************************
371     *                                   TOL                             *
372     *****************************************/
373     sciErr = getVarAddressFromPosition(pvApiCtx, 6, &piAddressVarSix);
374     if (sciErr.iErr)
375     {
376         printError(&sciErr, 0);
377         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 6);
378         return 0;
379     }
380
381     iErr = getScalarDouble(pvApiCtx, piAddressVarSix, &dblTOL);
382     if (iErr)
383     {
384         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be a real scalar.\n"), "eigs", 6, "opts.tol");
385         return 0;
386     }
387
388     if (C2F(isanan)(&dblTOL))
389     {
390         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be a real scalar.\n"), "eigs", 6, "opts.tol");
391         return 0;
392     }
393
394     /****************************************
395     *                                   NCV                             *
396     *****************************************/
397     sciErr = getVarAddressFromPosition(pvApiCtx, 7, &piAddressVarSeven);
398     if (sciErr.iErr)
399     {
400         printError(&sciErr, 0);
401         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 7);
402         return 0;
403     }
404
405     sciErr = getVarType(pvApiCtx, piAddressVarSeven, &TypeVarSeven);
406     if (sciErr.iErr || TypeVarSeven != sci_matrix)
407     {
408         printError(&sciErr, 0);
409         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv");
410         return 0;
411     }
412     else
413     {
414         if (isVarComplex(pvApiCtx, piAddressVarSeven))
415         {
416             Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv");
417             0;
418         }
419         else
420         {
421             sciErr = getVarDimension(pvApiCtx, piAddressVarSeven, &RowsSeven, &ColsSeven);
422             if (RowsSeven * ColsSeven > 1)
423             {
424                 Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv");
425                 return 0;
426             }
427
428             if (RowsSeven * ColsSeven == 1)
429             {
430                 sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarSeven, &RowsSeven, &ColsSeven, &dblNCV);
431                 if (sciErr.iErr)
432                 {
433                     printError(&sciErr, 0);
434                     Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 7);
435                     return 0;
436                 }
437
438                 if (dblNCV[0] != floor(dblNCV[0]))
439                 {
440                     Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv");
441                     return 0;
442                 }
443             }
444         }
445     }
446
447     /****************************************
448     *                           CHOLB                           *
449     *****************************************/
450     sciErr = getVarAddressFromPosition(pvApiCtx, 8, &piAddressVarEight);
451     if (sciErr.iErr)
452     {
453         printError(&sciErr, 0);
454         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 8);
455         return 0;
456     }
457
458     sciErr = getVarType(pvApiCtx, piAddressVarEight, &iTypeVarEight);
459     if (sciErr.iErr || iTypeVarEight != sci_matrix && iTypeVarEight != sci_boolean)
460     {
461         Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB");
462         return 0;
463     }
464
465     if (iTypeVarEight == sci_boolean)
466     {
467         iErr = getScalarBoolean(pvApiCtx, piAddressVarEight, &iCHOLB);
468         if (iErr)
469         {
470             Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB");
471             return 0;
472         }
473
474         if (iCHOLB != 1 && iCHOLB != 0)
475         {
476             Scierror(999, _("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB", "%f", "%t");
477             return 0;
478         }
479         dblCHOLB = (double) iCHOLB;
480     }
481
482     if (iTypeVarEight == sci_matrix)
483     {
484         iErr = getScalarDouble(pvApiCtx, piAddressVarEight, &dblCHOLB);
485         if (iErr)
486         {
487             Scierror(999, _("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB");
488             return 0;
489         }
490
491         if (dblCHOLB != 1 && dblCHOLB != 0)
492         {
493             Scierror(999, _("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB", "%f", "%t");
494             return 0;
495         }
496     }
497
498     if ( dblCHOLB ) // check that B is upper triangular with non zero element on the diagonal
499     {
500         if (!Bcomplex)
501         {
502             for (i = 0; i < N; i++)
503             {
504                 for (j = 0; j <= i; j++)
505                 {
506                     if (i == j && Breal[i + j * N] == 0)
507                     {
508                         Scierror(999, _("%s: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs");
509                         return 0;
510                     }
511                     else
512                     {
513                         if ( j < i && Breal[i + j * N] != 0 )
514                         {
515                             Scierror(999, _("%s: If opts.cholB is true, B should be upper triangular.\n"), "eigs");
516                             return 0;
517                         }
518                     }
519                 }
520             }
521         }
522         else
523         {
524             for (i = 0; i < N; i++)
525             {
526                 for (j = 0; j <= i; j++)
527                 {
528                     if (i == j && Bcplx[i + i * N].r == 0 && Bcplx[i + i * N].i == 0)
529                     {
530                         Scierror(999, _("%s: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs");
531                         return 0;
532                     }
533                     else
534                     {
535                         if ( j < i && (Bcplx[i + j * N].r != 0 || Bcplx[i + j * N].i != 0) )
536                         {
537                             Scierror(999, _("%s: If opts.cholB is true, B should be upper triangular.\n"), "eigs");
538                             return 0;
539                         }
540                     }
541                 }
542             }
543         }
544     }
545
546     /****************************************
547     *                           RESID                           *
548     *****************************************/
549     sciErr = getVarAddressFromPosition(pvApiCtx, 9, &piAddressVarNine);
550     if (sciErr.iErr)
551     {
552         printError(&sciErr, 0);
553         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 9);
554         return 0;
555     }
556
557     sciErr = getVarType(pvApiCtx, piAddressVarNine, &iTypeVarNine);
558     if (sciErr.iErr || iTypeVarNine != sci_matrix)
559     {
560         printError(&sciErr, 0);
561         Scierror(999, _("%s: Wrong type for input argument #%d: A real or complex matrix expected.\n"), "eigs", 9);
562         return 0;
563     }
564     else
565     {
566         sciErr = getVarDimension(pvApiCtx, piAddressVarNine, &iRowsNine, &iColsNine);
567         if (iRowsNine * iColsNine == 1 || iRowsNine * iColsNine != N)
568         {
569             Scierror(999, _("%s: Wrong dimension for input argument #%d: Start vector %s must be N by 1.\n"), "eigs", 9, "opts.resid");
570             return 0;
571         }
572     }
573
574     if (!Acomplex && !Bcomplex)
575     {
576         if (isVarComplex(pvApiCtx, piAddressVarNine))
577         {
578             Scierror(999, _("%s: Wrong type for input argument #%d: Start vector %s must be real for real problems.\n"), "eigs", 9, "opts.resid");
579             return 0;
580         }
581         else
582         {
583             sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarNine, &iRowsNine, &iColsNine, &RESID);
584             if (sciErr.iErr)
585             {
586                 printError(&sciErr, 0);
587                 Scierror(999, _("%s: Can not read input argument #%d.\n"), "eigs", 9);
588                 return 0;
589             }
590         }
591     }
592     else
593     {
594         sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddressVarNine, &iRowsNine, &iColsNine, &RESIDC);
595         if (sciErr.iErr)
596         {
597             printError(&sciErr, 0);
598             Scierror(999, _("%s: Can not read input argument #%d.\n"), "eigs", 9);
599             return 0;
600         }
601     }
602
603     /****************************************
604     *                           INFO                            *
605     *****************************************/
606     sciErr = getVarAddressFromPosition(pvApiCtx, 10, &piAddressVarTen);
607     if (sciErr.iErr)
608     {
609         printError(&sciErr, 0);
610         Scierror(999, _("%s: Can not read input argument #%d.\n"), "eigs", 9);
611         return 0;
612     }
613
614     iErr = getScalarInteger32(pvApiCtx, piAddressVarTen, &iINFO);
615     if (iErr)
616     {
617         Scierror(999, _("%s: Wrong type for input argument #%d: An integer expected.\n"), "eigs", 1);
618         return 0;
619     }
620
621     // Initialization output arguments
622     if (Lhs > 1)
623     {
624         RVEC = 1;
625     }
626
627     if (Acomplex || Bcomplex || !Asym)
628     {
629         eigenvalueC = (doublecomplex*)CALLOC((iNEV + 1), sizeof(doublecomplex));
630         if (RVEC)
631         {
632             eigenvectorC = (doublecomplex*)CALLOC(N * (iNEV + 1), sizeof(doublecomplex));
633         }
634     }
635     else
636     {
637         eigenvalue = (double*)CALLOC(iNEV, sizeof(double));
638         /* we should allocate eigenvector only if RVEC is true, but dseupd segfaults
639          if Z is not allocated even when RVEC is false, contrary to the docs.*/
640         eigenvector = (double*)CALLOC(iNEV * N, sizeof(double));
641     }
642
643     error = eigs(Areal, Acplx, N, Acomplex, Asym, Breal, Bcplx, Bcomplex, matB, iNEV, SIGMA, pstData, &dblMAXITER,
644                  &dblTOL, dblNCV, RESID, RESIDC, &iINFO, &dblCHOLB, INFO_EUPD, eigenvalue, eigenvector, eigenvalueC, eigenvectorC, RVEC);
645
646     switch (error)
647     {
648         case -1 :
649             if (Asym && !Acomplex && !Bcomplex)
650             {
651                 Scierror(999, _("%s: Wrong value for input argument #%d: For real symmetric problems, NCV must be k < NCV <= N.\n"), "eigs", 7);
652             }
653             else
654             {
655                 if (!Asym && !Acomplex && !Bcomplex)
656                 {
657                     Scierror(999, _("%s: Wrong value for input argument #%d: For real non symmetric problems, NCV must be k + 2 < NCV <= N.\n"), "eigs", 7);
658                 }
659                 else
660                 {
661                     Scierror(999, _("%s: Wrong value for input argument #%d: For complex problems, NCV must be k + 1 < NCV <= N.\n"), "eigs", 7);
662                 }
663             }
664             PutLhsVar();
665             return 0;
666
667         case -2 :
668             if (Asym && !Acomplex && !Bcomplex)
669             {
670                 Scierror(999, _("%s: Wrong value for input argument #%d: For real symmetric problems, k must be an integer in the range 1 to N - 1.\n"), "eigs", 3);
671             }
672             else
673             {
674                 Scierror(999, _("%s: Wrong value for input argument #%d: For real non symmetric or complex problems, k must be an integer in the range 1 to N - 2.\n"), "eigs", 3);
675             }
676             PutLhsVar();
677             return 0;
678
679         case -3 :
680             Scierror(999, _("%s: Error with input argument #%d: B is not positive definite. Try with sigma='SM' or sigma=scalar.\n"), "eigs", 2);
681             ReturnArguments(pvApiCtx);
682             return 0;
683
684         case -4 :
685             if (!Acomplex && !Bcomplex)
686             {
687                 if (Asym)
688                 {
689                     Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "DSAUPD", iINFO);
690                 }
691                 else
692                 {
693                     Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "DNAUPD", iINFO);
694                 }
695             }
696             else
697             {
698                 Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "ZNAUPD", iINFO);
699             }
700             PutLhsVar();
701             return 0;
702
703         case -5 :
704             if (!Acomplex && !Bcomplex)
705             {
706                 if (Asym)
707                 {
708                     Scierror(999, _("%s: Error with %s: unknown mode returned.\n"), "eigs", "DSAUPD");
709                 }
710                 else
711                 {
712                     Scierror(999, _("%s: Error with %s: unknown mode returned.\n"), "eigs", "DNAUPD");
713                 }
714             }
715             else
716             {
717                 Scierror(999, _("%s: Error with %s: unknown mode returned.\n"), "eigs", "ZNAUPD");
718             }
719             PutLhsVar();
720             return 0;
721
722         case -6 :
723             if (!Acomplex && !Bcomplex)
724             {
725                 if (Asym)
726                 {
727                     Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "DSEUPD", INFO_EUPD);
728                 }
729                 else
730                 {
731                     Scierror(999, _("%s: Error with %s: info = %d \n"), "eigs", "DNEUPD", INFO_EUPD);
732                 }
733             }
734             else
735             {
736                 Scierror(999,  _("%s: Error with %s: info = %d \n"), "eigs", "ZNEUPD", INFO_EUPD);
737             }
738             PutLhsVar();
739             return 0;
740
741         case -7 :
742             Scierror(999, _("%s: A - sigma * B is not inversible, try with a different value of sigma.\n"), "eigs");
743             PutLhsVar();
744             return 0;
745     }
746
747     if (Lhs <= 1)
748     {
749         if (eigenvalue)
750         {
751             sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, 1, eigenvalue);
752             FREE(eigenvalue);
753             FREE(eigenvector);
754         }
755         else if (eigenvalueC)
756         {
757             sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, 1, eigenvalueC);
758             FREE(eigenvalueC);
759         }
760         if (sciErr.iErr)
761         {
762             printError(&sciErr, 0);
763             Scierror(999, _("%s: Memory allocation error.\n"), fname);
764             return 0;
765         }
766
767         LhsVar(1) = Rhs + 1;
768     }
769     else
770     {
771         // create a matrix which contains the eigenvalues
772         if (eigenvalue)
773         {
774             mat_eigenvalue = (double*)CALLOC(iNEV * iNEV, sizeof(double));
775             for (i = 0; i < iNEV; i++)
776             {
777                 mat_eigenvalue[i * iNEV + i] = eigenvalue[i];
778             }
779             sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, iNEV, mat_eigenvalue);
780             FREE(eigenvalue);
781             FREE(mat_eigenvalue);
782         }
783         else if (eigenvalueC)
784         {
785             mat_eigenvalueC = (doublecomplex*)CALLOC(iNEV * iNEV, sizeof(doublecomplex));
786             for (i = 0; i < iNEV; i++)
787             {
788                 mat_eigenvalueC[i * iNEV + i] = eigenvalueC[i];
789             }
790             sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 1, iNEV, iNEV, mat_eigenvalueC);
791             FREE(eigenvalueC);
792             FREE(mat_eigenvalueC);
793         }
794
795         if (sciErr.iErr)
796         {
797             printError(&sciErr, 0);
798             Scierror(999, _("%s: Memory allocation error.\n"), fname);
799             return 0;
800         }
801
802         if (eigenvector)
803         {
804             sciErr = createMatrixOfDouble(pvApiCtx, Rhs + 2, N, iNEV, eigenvector);
805             FREE(eigenvector);
806         }
807         else if (eigenvectorC)
808         {
809             sciErr = createComplexZMatrixOfDouble(pvApiCtx, Rhs + 2, N, iNEV, eigenvectorC);
810             FREE(eigenvectorC);
811         }
812         if (sciErr.iErr)
813         {
814             printError(&sciErr, 0);
815             Scierror(999, _("%s: Memory allocation error.\n"), fname);
816             return 0;
817         }
818
819         LhsVar(1) = Rhs + 1;
820         LhsVar(2) = Rhs + 2;
821     }
822
823     if (iTypeVarFour == sci_strings)
824     {
825         freeAllocatedSingleString(pstData);
826     }
827
828     if (matB != 0)
829     {
830         if (Acomplex && !Bcomplex)
831         {
832             free(Bcplx);
833         }
834         if (!Acomplex && Bcomplex)
835         {
836             free(Acplx);
837         }
838     }
839
840     PutLhsVar();
841     return 0;
842 }