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