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