Merge remote-tracking branch 'origin/master' into YaSp
[scilab.git] / scilab / modules / arnoldi / sci_gateway / c / sci_zneupd.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) ????-2008 - INRIA
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 "api_scilab.h"
16 #include "gw_arnoldi.h"
17 #include "core_math.h"
18 #include "doublecomplex.h"
19 #include "localization.h"
20 #include "Scierror.h"
21 /*--------------------------------------------------------------------------*/
22 extern int C2F(zneupd)(int * rvec, char * howmny, int * select,
23                        doublecomplex * d, doublecomplex * z, int * ldz,
24                        doublecomplex * sigma, doublecomplex * workev,
25                        char * bmat, int * n, char * which, int * nev,
26                        double *  tol, doublecomplex * resid, int * ncv,
27                        doublecomplex * v, int * ldv, int * iparam, int * ipntr,
28                        doublecomplex * workd, doublecomplex * workl,
29                        int * lworkl, double * rwork, int * info);
30 /*--------------------------------------------------------------------------*/
31 int sci_zneupd(char *fname, void *pvApiCtx)
32 {
33     SciErr sciErr;
34
35     int* piAddrpRVEC    = NULL;
36     int* pRVEC          = NULL;
37     int* piAddrpHOWMANY = NULL;
38     char* pHOWMANY      = NULL;
39     int* piAddrpSELECT  = NULL;
40     int* pSELECT        = NULL;
41     int* piAddrpBMAT    = NULL;
42     char* pBMAT         = NULL;
43     int* piAddrpN       = NULL;
44     int* pN             = NULL;
45     int* piAddrpWHICH   = NULL;
46     char* pWHICH        = NULL;
47     int* piAddrpNEV     = NULL;
48     int* pNEV           = NULL;
49     int* piAddrpTOL     = NULL;
50     double* pTOL        = NULL;
51     int* piAddrpNCV     = NULL;
52     int* pNCV           = NULL;
53     int* piAddrpIPARAM  = NULL;
54     int* pIPARAM        = NULL;
55     int* piAddrpIPNTR   = NULL;
56     int* pIPNTR         = NULL;
57     int* piAddrpRWORK   = NULL;
58     double* pRWORK      = NULL;
59     int* piAddrpINFO    = NULL;
60     int* pINFO          = NULL;
61
62     int* piAddrpD           = NULL;
63     doublecomplex* pD       = NULL;
64     int* piAddrpZ           = NULL;
65     doublecomplex* pZ       = NULL;
66     int* piAddrpSIGMA       = NULL;
67     doublecomplex* pSIGMA   = NULL;
68     int* piAddrpWORKev      = NULL;
69     doublecomplex* pWORKev  = NULL;
70     int* piAddrpRESID       = NULL;
71     doublecomplex* pRESID   = NULL;
72     int* piAddrpWORKD       = NULL;
73     doublecomplex* pV       = NULL;
74     int* piAddrpV           = NULL;
75     doublecomplex* pWORKD   = NULL;
76     int* piAddrpWORKL       = NULL;
77     doublecomplex* pWORKL   = NULL;
78
79     int mRVEC,     nRVEC;
80     int mHOWMANY,  nHOWMANY;
81     int mSELECT,   nSELECT;
82     int D,        mD,        nD;
83     int Z,        mZ,        nZ;
84     int mSIGMA,    nSIGMA;
85     int mWORKev,   nWORKev;
86     int mBMAT,     nBMAT;
87     int mN,        nN;
88     int mWHICH,    nWHICH;
89     int mNEV,      nNEV;
90     int mTOL,      nTOL;
91     int RESID,    mRESID,    nRESID;
92     int mNCV,      nNCV;
93     int mV,        nV;
94     int IPARAM,   mIPARAM,   nIPARAM;
95     int IPNTR,    mIPNTR,    nIPNTR;
96     int WORKD,    mWORKD,    nWORKD;
97     int WORKL,    mWORKL,    nWORKL;
98     int RWORK,    mRWORK,    nRWORK;
99     int INFO,     mINFO,     nINFO;
100
101     int minlhs = 1, minrhs = 21, maxlhs = 9, maxrhs = 21;
102     int LDZ, LDV, LWORKL;
103     int sizeWORKL = 0;
104
105     CheckInputArgument(pvApiCtx, minrhs, maxrhs);
106     CheckOutputArgument(pvApiCtx, minlhs, maxlhs);
107
108     /*                                                  VARIABLE = NUMBER   */
109     sciErr = getVarAddressFromPosition(pvApiCtx, 1, &piAddrpRVEC);
110     if (sciErr.iErr)
111     {
112         printError(&sciErr, 0);
113         return 1;
114     }
115
116     // Retrieve a matrix of double at position 1.
117     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpRVEC, &mRVEC, &nRVEC, &pRVEC);
118     if (sciErr.iErr)
119     {
120         printError(&sciErr, 0);
121         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 1);
122         return 1;
123     }
124
125     sciErr = getVarAddressFromPosition(pvApiCtx, 3, &piAddrpSELECT);
126     if (sciErr.iErr)
127     {
128         printError(&sciErr, 0);
129         return 1;
130     }
131
132     // Retrieve a matrix of double at position 3.
133     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpSELECT, &mSELECT, &nSELECT, &pSELECT);
134     if (sciErr.iErr)
135     {
136         printError(&sciErr, 0);
137         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 3);
138         return 1;
139     }
140
141     sciErr = getVarAddressFromPosition(pvApiCtx, 4, &piAddrpD);
142     if (sciErr.iErr)
143     {
144         printError(&sciErr, 0);
145         return 1;
146     }
147
148     // Retrieve a matrix of double at position 4.
149     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpD, &mD, &nD, &pD);
150     if (sciErr.iErr)
151     {
152         printError(&sciErr, 0);
153         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 4);
154         return 1;
155     }
156     D =  4;
157
158     sciErr = getVarAddressFromPosition(pvApiCtx, 5, &piAddrpZ);
159     if (sciErr.iErr)
160     {
161         printError(&sciErr, 0);
162         return 1;
163     }
164
165     // Retrieve a matrix of double at position 5.
166     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpZ, &mZ, &nZ, &pZ);
167     if (sciErr.iErr)
168     {
169         printError(&sciErr, 0);
170         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 5);
171         return 1;
172     }
173     Z =  5;
174
175     sciErr = getVarAddressFromPosition(pvApiCtx, 6, &piAddrpSIGMA);
176     if (sciErr.iErr)
177     {
178         printError(&sciErr, 0);
179         return 1;
180     }
181
182     // Retrieve a matrix of double at position 6.
183     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpSIGMA, &mSIGMA, &nSIGMA, &pSIGMA);
184     if (sciErr.iErr)
185     {
186         printError(&sciErr, 0);
187         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 6);
188         return 1;
189     }
190
191     sciErr = getVarAddressFromPosition(pvApiCtx, 7, &piAddrpWORKev);
192     if (sciErr.iErr)
193     {
194         printError(&sciErr, 0);
195         return 1;
196     }
197
198     // Retrieve a matrix of double at position 7.
199     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpWORKev, &mWORKev, &nWORKev, &pWORKev);
200     if (sciErr.iErr)
201     {
202         printError(&sciErr, 0);
203         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 7);
204         return 1;
205     }
206
207     sciErr = getVarAddressFromPosition(pvApiCtx, 9, &piAddrpN);
208     if (sciErr.iErr)
209     {
210         printError(&sciErr, 0);
211         return 1;
212     }
213
214     // Retrieve a matrix of double at position 9.
215     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpN, &mN, &nN, &pN);
216     if (sciErr.iErr)
217     {
218         printError(&sciErr, 0);
219         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 9);
220         return 1;
221     }
222
223     sciErr = getVarAddressFromPosition(pvApiCtx, 11, &piAddrpNEV);
224     if (sciErr.iErr)
225     {
226         printError(&sciErr, 0);
227         return 1;
228     }
229
230     // Retrieve a matrix of double at position 11.
231     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpNEV, &mNEV, &nNEV, &pNEV);
232     if (sciErr.iErr)
233     {
234         printError(&sciErr, 0);
235         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 11);
236         return 1;
237     }
238
239     sciErr = getVarAddressFromPosition(pvApiCtx, 12, &piAddrpTOL);
240     if (sciErr.iErr)
241     {
242         printError(&sciErr, 0);
243         return 1;
244     }
245
246     // Retrieve a matrix of double at position 12.
247     sciErr = getMatrixOfDouble(pvApiCtx, piAddrpTOL, &mTOL, &nTOL, &pTOL);
248     if (sciErr.iErr)
249     {
250         printError(&sciErr, 0);
251         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 12);
252         return 1;
253     }
254
255     sciErr = getVarAddressFromPosition(pvApiCtx, 13, &piAddrpRESID);
256     if (sciErr.iErr)
257     {
258         printError(&sciErr, 0);
259         return 1;
260     }
261
262     // Retrieve a matrix of double at position 13.
263     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpRESID, &mRESID, &nRESID, &pRESID);
264     if (sciErr.iErr)
265     {
266         printError(&sciErr, 0);
267         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 13);
268         return 1;
269     }
270     RESID = 13;
271
272     sciErr = getVarAddressFromPosition(pvApiCtx, 14, &piAddrpNCV);
273     if (sciErr.iErr)
274     {
275         printError(&sciErr, 0);
276         return 1;
277     }
278
279     // Retrieve a matrix of double at position 14.
280     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpNCV, &mNCV, &nNCV, &pNCV);
281     if (sciErr.iErr)
282     {
283         printError(&sciErr, 0);
284         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 14);
285         return 1;
286     }
287
288     sciErr = getVarAddressFromPosition(pvApiCtx, 15, &piAddrpV);
289     if (sciErr.iErr)
290     {
291         printError(&sciErr, 0);
292         return 1;
293     }
294
295     // Retrieve a matrix of double at position 15.
296     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpV, &mV, &nV, &pV);
297     if (sciErr.iErr)
298     {
299         printError(&sciErr, 0);
300         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 15);
301         return 1;
302     }
303
304     sciErr = getVarAddressFromPosition(pvApiCtx, 16, &piAddrpIPARAM);
305     if (sciErr.iErr)
306     {
307         printError(&sciErr, 0);
308         return 1;
309     }
310
311     // Retrieve a matrix of double at position 16.
312     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpIPARAM, &mIPARAM, &nIPARAM, &pIPARAM);
313     if (sciErr.iErr)
314     {
315         printError(&sciErr, 0);
316         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 16);
317         return 1;
318     }
319
320     IPARAM = 16;
321     sciErr = getVarAddressFromPosition(pvApiCtx, 17, &piAddrpIPNTR);
322     if (sciErr.iErr)
323     {
324         printError(&sciErr, 0);
325         return 1;
326     }
327
328     // Retrieve a matrix of double at position 17.
329     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpIPNTR, &mIPNTR, &nIPNTR, &pIPNTR);
330     if (sciErr.iErr)
331     {
332         printError(&sciErr, 0);
333         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 17);
334         return 1;
335     }
336
337     IPNTR = 17;
338
339     sciErr = getVarAddressFromPosition(pvApiCtx, 18, &piAddrpWORKD);
340     if (sciErr.iErr)
341     {
342         printError(&sciErr, 0);
343         return 1;
344     }
345
346     // Retrieve a matrix of double at position 18.
347     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpWORKD, &mWORKD, &nWORKD, &pWORKD);
348     if (sciErr.iErr)
349     {
350         printError(&sciErr, 0);
351         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 18);
352         return 1;
353     }
354     WORKD = 18;
355
356     sciErr = getVarAddressFromPosition(pvApiCtx, 19, &piAddrpWORKL);
357     if (sciErr.iErr)
358     {
359         printError(&sciErr, 0);
360         return 1;
361     }
362
363     // Retrieve a matrix of double at position 19.
364     sciErr = getComplexZMatrixOfDouble(pvApiCtx, piAddrpWORKL, &mWORKL, &nWORKL, &pWORKL);
365     if (sciErr.iErr)
366     {
367         printError(&sciErr, 0);
368         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 19);
369         return 1;
370     }
371     WORKL = 19;
372
373     sciErr = getVarAddressFromPosition(pvApiCtx, 20, &piAddrpRWORK);
374     if (sciErr.iErr)
375     {
376         printError(&sciErr, 0);
377         return 1;
378     }
379
380     // Retrieve a matrix of double at position 20.
381     sciErr = getMatrixOfDouble(pvApiCtx, piAddrpRWORK, &mRWORK, &nRWORK, &pRWORK);
382     if (sciErr.iErr)
383     {
384         printError(&sciErr, 0);
385         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 20);
386         return 1;
387     }
388
389     RWORK = 20;
390     sciErr = getVarAddressFromPosition(pvApiCtx, 21, &piAddrpINFO);
391     if (sciErr.iErr)
392     {
393         printError(&sciErr, 0);
394         return 1;
395     }
396
397     // Retrieve a matrix of double at position 21.
398     sciErr = getMatrixOfDoubleAsInteger(pvApiCtx, piAddrpINFO, &mINFO, &nINFO, &pINFO);
399     if (sciErr.iErr)
400     {
401         printError(&sciErr, 0);
402         Scierror(202, _("%s: Wrong type for argument %d: A real expected.\n"), fname, 21);
403         return 1;
404     }
405
406     INFO = 21;
407
408     LWORKL = mWORKL * nWORKL;
409     LDV = Max(1, pN[0]);
410     LDZ = LDV;
411
412     /* Check some sizes */
413     if (mIPARAM*nIPARAM != 11)
414     {
415         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "IPARAM", 11);
416         return 1;
417     }
418
419     if (mIPNTR*nIPNTR != 14)
420     {
421         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "IPNTR", 14);
422         return 1;
423     }
424
425     if (mRESID*nRESID != pN[0])
426     {
427         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "RESID", pN[0]);
428         return 1;
429     }
430
431     if (mWORKD * nWORKD < 3 * pN[0])
432     {
433         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "WORKD", 3 * pN[0]);
434         return 1;
435     }
436
437     if (mSELECT*nSELECT != pNCV[0])
438     {
439         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "SELECT", pNCV[0]);
440         return 1;
441     }
442
443     if (mD*nD != (pNEV[0] + 1))
444     {
445         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "D", pNEV[0] + 1);
446         return 1;
447     }
448
449     if ((mZ != pN[0]) || (nZ != pNEV[0]))
450     {
451         Scierror(999, _("%s: Wrong size for input argument %s: A matrix of size %dx%d expected.\n"), fname, "Z", pN[0], pNEV[0]);
452         return 1;
453     }
454
455     if (mWORKev*nWORKev != 2 * pNCV[0])
456     {
457         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "WORKev", 2 * pNCV[0]);
458         return 1;
459     }
460
461     if ((mV != pN[0]) || (nV != pNCV[0]))
462     {
463         Scierror(999, _("%s: Wrong size for input argument %s: A matrix of size %dx%d expected.\n"), fname, "V", pN[0], pNCV[0]);
464         return 1;
465     }
466
467     sizeWORKL = 3 * pNCV[0] * pNCV[0] + 5 * pNCV[0];
468
469     if ((mWORKL * nWORKL < sizeWORKL))
470     {
471         Scierror(999, _("%s: Wrong size for input argument %s: An array of size %d expected.\n"), fname, "WORKL", sizeWORKL);
472         return 1;
473     }
474
475     sciErr = getVarAddressFromPosition(pvApiCtx, 2, &piAddrpHOWMANY);
476     if (sciErr.iErr)
477     {
478         printError(&sciErr, 0);
479         return 1;
480     }
481
482     // Retrieve a matrix of double at position 2.
483     if (getAllocatedSingleString(pvApiCtx, piAddrpHOWMANY, &pHOWMANY))
484     {
485         Scierror(202, _("%s: Wrong type for argument #%d: A string expected.\n"), fname, 2);
486         return 1;
487     }
488
489     sciErr = getVarAddressFromPosition(pvApiCtx, 8, &piAddrpBMAT);
490     if (sciErr.iErr)
491     {
492         freeAllocatedSingleString(pHOWMANY);
493         printError(&sciErr, 0);
494         return 1;
495     }
496
497     // Retrieve a matrix of double at position 8.
498     if (getAllocatedSingleString(pvApiCtx, piAddrpBMAT, &pBMAT))
499     {
500         freeAllocatedSingleString(pHOWMANY);
501         Scierror(202, _("%s: Wrong type for argument #%d: A string expected.\n"), fname, 8);
502         return 1;
503     }
504
505     sciErr = getVarAddressFromPosition(pvApiCtx, 10, &piAddrpWHICH);
506     if (sciErr.iErr)
507     {
508         freeAllocatedSingleString(pHOWMANY);
509         freeAllocatedSingleString(pBMAT);
510         printError(&sciErr, 0);
511         return 1;
512     }
513
514     // Retrieve a matrix of double at position 10.
515     if (getAllocatedSingleString(pvApiCtx, piAddrpWHICH, &pWHICH))
516     {
517         freeAllocatedSingleString(pHOWMANY);
518         freeAllocatedSingleString(pBMAT);
519         Scierror(202, _("%s: Wrong type for argument #%d: A string expected.\n"), fname, 10);
520         return 1;
521     }
522
523     C2F(zneupd)(pRVEC, pHOWMANY, pSELECT, pD,
524                 pZ, &LDZ, pSIGMA, pWORKev,
525                 pBMAT, pN, pWHICH, pNEV,
526                 pTOL, pRESID, pNCV, pV,
527                 &LDV, pIPARAM, pIPNTR, pWORKD,
528                 pWORKL, &LWORKL, pRWORK, pINFO);
529
530     freeAllocatedSingleString(pHOWMANY);
531     freeAllocatedSingleString(pBMAT);
532     freeAllocatedSingleString(pWHICH);
533
534     if (pINFO[0] < 0)
535     {
536         C2F(errorinfo)("zneupd", pINFO, 6L);
537         return 0;
538     }
539
540     AssignOutputVariable(pvApiCtx, 1) = D;
541     AssignOutputVariable(pvApiCtx, 2) = Z;
542     AssignOutputVariable(pvApiCtx, 3) = RESID;
543     AssignOutputVariable(pvApiCtx, 4) = IPARAM;
544     AssignOutputVariable(pvApiCtx, 5) = IPNTR;
545     AssignOutputVariable(pvApiCtx, 6) = WORKD;
546     AssignOutputVariable(pvApiCtx, 7) = WORKL;
547     AssignOutputVariable(pvApiCtx, 8) = RWORK;
548     AssignOutputVariable(pvApiCtx, 9) = INFO;
549
550     ReturnArguments(pvApiCtx);
551
552     return 0;
553 }
554 /*--------------------------------------------------------------------------*/