1 #include "mex.h"
2 #include "stack-c.h"
3 #include "sci_gateway.h"
4 #include <stdlib.h>
5 #include <stdio.h>
6 #include <string.h>
7
8 #define MUMPS_ARITH_d 2
9 #define MUMPS_ARITH_z 8
10
11 #if MUMPS_ARITH == MUMPS_ARITH_z
12
13 # include "zmumps_c.h"
14 # define dmumps_c zmumps_c
15 # define dmumps_par zmumps_par
16 # define DMUMPS_STRUC_C ZMUMPS_STRUC_C
17 # define DMUMPS_alloc ZMUMPS_alloc
18 # define DMUMPS_free ZMUMPS_free
19 # define double2 mumps_double_complex
20
21 #elif MUMPS_ARITH == MUMPS_ARITH_d
22
23 # include "dmumps_c.h"
24 # define double2 double
25 # define EXTRACT_CMPLX_FROM_C_TO_SCILAB EXTRACT_DOUBLE_FROM_C_TO_SCILAB
26 # define EXTRACT_CMPLX_FROM_SCILAB_TOPTR EXTRACT_FROM_SCILAB_TOPTR
27
28 #else
29
30 # error "Only d and z arithmetics are supported"
31
32 #endif
33
34
35 #define nb_RHS 12
36
37 #define MYFREE(ptr)\
38 if(ptr){ \
39 free(ptr); \
40 ptr=0;} \
41
42
43 #define EXTRACT_FROM_SCILAB_TOPTR(it,ptr_scilab1,ptr_scilab2,mumpspointer,type,length)\
44 if(ptr_scilab1[0] != -9999){ \
45 free(mumpspointer); \
46 mumpspointer = (type *) malloc(length*sizeof(type)); \
47 for(i=0;i<length;i++){ \
48 mumpspointer[i] = ptr_scilab1[i]; \
49 } \
50 }
51
52 #define EXTRACT_FROM_SCILAB_TOARR(ptr_scilab,mumpsarray,type,length) \
53 if(ptr_scilab[0] != -9999){ \
54 for(i=0;i<length;i++){ \
55 if(ptr_scilab[i] != -9998){ \
56 mumpsarray[i] = ptr_scilab[i]; \
57 } \
58 } \
59 }
60
61 #define EXTRACT_INT_FROM_C_TO_SCILAB(num,ptr_scilab,mumpsptr,length1,length2,one)\
62 if(mumpsptr == 0){ \
63 CreateVar(nb_RHS+num,"i",&one,&one,&ptr_scilab); \
64 *istk(ptr_scilab)=-9999; \
65 }else{ \
66 CreateVar(nb_RHS+num,"i",&length1,&length2,&ptr_scilab); \
67 int l=length1*length2; \
68 for (i=0;i<l;i++){ \
69 *istk(ptr_scilab+i)=(mumpsptr)[i];} \
70 } \
71 LhsVar(num)=nb_RHS+num;
72
73 #define EXTRACT_DOUBLE_FROM_C_TO_SCILAB(num,it,ptr_scilab1,ptr_scilab2,mumpsptr,length1,length2,one)\
74 if(mumpsptr == 0){ \
75 CreateVar(nb_RHS+num,"d",&one,&one,&ptr_scilab1); \
76 *stk(ptr_scilab1)=-9999; \
77 }else{ \
78 CreateVar(nb_RHS+num,"d",&length1,&length2,&ptr_scilab1); \
79 int l=length1*length2; \
80 for (i=0;i<l;i++){ \
81 *stk(ptr_scilab1+i)=(mumpsptr)[i];} \
82 } \
83 LhsVar(num)=nb_RHS+num;
84
85 #if MUMPS_ARITH == MUMPS_ARITH_z
86
87 #define EXTRACT_CMPLX_FROM_SCILAB_TOPTR(it,ptr_scilab1,ptr_scilab2,mumpsptr,type,length)\
88 if(ptr_scilab1[0] != -9999){ \
89 free(mumpsptr); \
90 mumpsptr=(double2 *) malloc(length*sizeof(double2)); \
91 for(i=0;i<length;i++){ \
92 (mumpsptr[i]).r = ptr_scilab1[i];} \
93 if(it==1){ \
94 for(i=0;i<length;i++){ \
95 (mumpsptr[i]).i = ptr_scilab2[i];} \
96 }else{ \
97 for(i=0;i<length;i++){ \
98 (mumpsptr[i]).i = 0.0;} \
99 } \
100 }
101
102 #define EXTRACT_CMPLX_FROM_C_TO_SCILAB(num,it,ptr_scilab1,ptr_scilab2,mumpsptr,length1,length2,one)\
103 if(it!=1){ \
104 Scierror(999,"Internal error, the variable built must be complex.\n");} \
105 if(mumpsptr == 0){ \
106 CreateCVar(nb_RHS+num,"d",&it, &one,&one,&ptr_scilab1,&ptr_scilab2); \
107 *stk(ptr_scilab1) = -9999; \
108 *stk(ptr_scilab2) = -9999; \
109 }else{ \
110 int l=length1*length2; \
111 CreateCVar(nb_RHS+num,"d",&it,&length1,&length2,&ptr_scilab1,&ptr_scilab2); \
112 for(i=0;i<l;i++){ \
113 *stk(ptr_scilab1+i) = (mumpsptr[i]).r; \
114 } \
115 for(i=0;i<l;i++){ \
116 *stk(ptr_scilab2+i) = (mumpsptr[i]).i; \
117 } \
118 } \
119 LhsVar(num)=nb_RHS+num;
120
121 #endif
122
123
124
DMUMPS_free(DMUMPS_STRUC_C ** dmumps_par)125 void DMUMPS_free(DMUMPS_STRUC_C **dmumps_par){
126 if(*dmumps_par){
127 MYFREE( (*dmumps_par)->irn );
128 MYFREE( (*dmumps_par)->jcn );
129 MYFREE( (*dmumps_par)->a );
130 MYFREE( (*dmumps_par)->irn_loc );
131 MYFREE( (*dmumps_par)->jcn_loc );
132 MYFREE( (*dmumps_par)->a_loc );
133 MYFREE( (*dmumps_par)->eltptr );
134 MYFREE( (*dmumps_par)->eltvar );
135 MYFREE( (*dmumps_par)->a_elt );
136 MYFREE( (*dmumps_par)->perm_in );
137 MYFREE( (*dmumps_par)->colsca );
138 MYFREE( (*dmumps_par)->rowsca );
139 MYFREE( (*dmumps_par)->pivnul_list );
140 MYFREE( (*dmumps_par)->listvar_schur );
141 MYFREE( (*dmumps_par)->sym_perm );
142 MYFREE( (*dmumps_par)->uns_perm );
143 MYFREE( (*dmumps_par)->irhs_ptr);
144 MYFREE( (*dmumps_par)->irhs_sparse);
145 MYFREE( (*dmumps_par)->rhs_sparse);
146 MYFREE( (*dmumps_par)->rhs);
147 MYFREE( (*dmumps_par)->redrhs);
148 MYFREE(*dmumps_par);
149 }
150 }
151
DMUMPS_alloc(DMUMPS_STRUC_C ** dmumps_par)152 void DMUMPS_alloc(DMUMPS_STRUC_C **dmumps_par){
153
154 *dmumps_par = (DMUMPS_STRUC_C *) malloc(sizeof(DMUMPS_STRUC_C));
155 (*dmumps_par)->irn = NULL;
156 (*dmumps_par)->jcn = NULL;
157 (*dmumps_par)->a = NULL;
158 (*dmumps_par)->irn_loc = NULL;
159 (*dmumps_par)->jcn_loc = NULL;
160 (*dmumps_par)->a_loc = NULL;
161 (*dmumps_par)->eltptr = NULL;
162 (*dmumps_par)->eltvar = NULL;
163 (*dmumps_par)->a_elt = NULL;
164 (*dmumps_par)->perm_in = NULL;
165 (*dmumps_par)->colsca = NULL;
166 (*dmumps_par)->rowsca = NULL;
167 (*dmumps_par)->rhs = NULL;
168 (*dmumps_par)->redrhs = NULL;
169 (*dmumps_par)->irhs_ptr = NULL;
170 (*dmumps_par)->irhs_sparse = NULL;
171 (*dmumps_par)->rhs_sparse = NULL;
172 (*dmumps_par)->pivnul_list = NULL;
173 (*dmumps_par)->listvar_schur = NULL;
174 (*dmumps_par)->schur = NULL;
175 (*dmumps_par)->sym_perm = NULL;
176 (*dmumps_par)->uns_perm = NULL;
177 }
178
179
180
dmumpsc(char * fname)181 static int dmumpsc(char *fname){
182
183
184 /* RhsVar parameters */
185 int njob, mjob, ljob, mint, nint, lint, nsym, msym, lsym, nA, mA, nRHS, nREDRHS, mRHS,lRHS, liRHS;
186 int mREDRHS,lREDRHS,liREDRHS;
187 int nicntl, micntl, licntl, ncntl, mcntl, lcntl, nperm, mperm, lperm;
188 int ncols, mcols, lcols, licols, nrows, mrows, lrows, lirows, ns_schu , ms_schu, ls_schu;
189 int nv_schu, mv_schu, lv_schu, nschu, mschu, lschu;
190 int type_rhs, mtype_rhs, ntype_rhs, ltype_rhs;
191
192 /* LhsVar parameters */
193 int linfog, lrinfog, lrhsout,lrhsouti, linstout, lschurout, lschurouti, ldef;
194 int lpivnul_list, lmapp, lsymperm, lunsperm;
195 int one=1, temp1=40, temp2=40, temp3, temp4;
196 int it, itRHS, itREDRHS; /* parameter for real/complex types */
197
198 int i,j,k1,k2, nb_in_row,netrue;
199 int *ptr_int;
200 double *ptr_double;
201 double *ptr_scilab;
202 #if MUMPS_ARITH == MUMPS_ARITH_z
203 double * ptri_scilab;
204 #endif
205
206 /* Temporary length variables */
207 int len1, len2;
208 /* Temporary pointers in stack */
209 int stkptr, stkptri;
210
211 /* C pointer for input parameters */
212 int inst_address;
213 int ne,inst;
214 int *irn_in,*jcn_in;
215
216 /* Variable for multiple and sparse RHS*/
217 int posrhs, posschur, nz_RHS,col_ind,k;
218 int *irhs_ptr;
219 int *irhs_sparse;
220 double *rhs_sparse;
221 #if MUMPS_ARITH == MUMPS_ARITH_z
222 double *im_rhs_sparse;
223 char * function_name="zmumpsc";
224 #else
225 char * function_name="dmumpsc";
226 #endif
227
228 SciSparse A;
229 SciSparse RHS_SPARSE;
230 DMUMPS_STRUC_C *dmumps_par;
231
232 int dosolve=0;
233 int donullspace=0;
234 int doanal = 0;
235 /* Check number of input parameters */
236 CheckRhs(11,12);
237
238 /* Get job value. njob/mjob are the dimensions of variable job. */
239 GetRhsVar(2,"i",&mjob,&njob,&ljob);
240 dosolve = (*istk(ljob) == 3 || *istk(ljob) == 5 ||*istk(ljob) == 6);
241 doanal = (*istk(ljob) == 1 || *istk(ljob) == 4 || *istk(ljob) == 6);
242 if(*istk(ljob) == -1){
243
244 DMUMPS_alloc(&dmumps_par);
245 GetRhsVar(1,"i",&msym,&nsym,&lsym);
246 dmumps_par->sym=*istk(lsym);
247 dmumps_par->job = -1;
248 dmumps_par->par = 1;
249 dmumps_c(dmumps_par);
250 dmumps_par->nz = -1;
251 dmumps_par->nz_alloc=-1;
252 it=1;
253 }else{
254 /* Obtain pointer on instance */
255 GetRhsVar(10,"i",&mint,&nint,&lint);
256 inst_address=*istk(lint); /* EXTRACT_FROM_SCILAB_TOVAL(INST,inst_address); */
257 ptr_int = (int *) inst_address;
258
259 dmumps_par = (DMUMPS_STRUC_C *) ptr_int;
260 if(*istk(ljob) == -2){
261 dmumps_par->job = -2;
262 dmumps_c(dmumps_par);
263 DMUMPS_free(&dmumps_par);
264 }else{
265 /* Get the sparse matrix A */
266 GetRhsVar(12,"s",&mA,&nA,&A);
267
268 if (nA != mA || mA<1 ){
269 Scierror(999,"%s: Bad dimensions for mat\n",function_name);
270 return 0;
271 }
272
273 ne=A.nel;
274 dmumps_par->n = nA;
275
276 if(dmumps_par->sym != 0){
277 netrue = (nA+ne)/2;
278 }else{
279 netrue = ne;
280 }
281
282 if(dmumps_par->nz_alloc < netrue ||dmumps_par->nz_alloc >= 2*netrue){
283 MYFREE(dmumps_par->jcn);
284 MYFREE(dmumps_par->irn);
285 MYFREE(dmumps_par->a);
286
287 dmumps_par->jcn = (int*)malloc(netrue*sizeof(int));
288 dmumps_par->irn = (int*)malloc(netrue*sizeof(int));
289 dmumps_par->a = (double2 *) malloc(netrue*sizeof(double2));
290 dmumps_par->nz_alloc = netrue;
291 }
292 /* Check for symmetry in order to initialize only
293 * lower triangle on entry to symmetric MUMPS code */
294 if ((dmumps_par->sym)==0){
295 /*
296 * Unsymmetric case:
297 * build irn from mnel for MUMPS format
298 * mA : number of rows
299 */
300
301 if(doanal){
302 for(i=0;i<ne;i++){
303 (dmumps_par->jcn)[i]=(A.icol)[i];}
304 k1=0;
305 for (k2=1;k2<mA+1;k2++){
306 nb_in_row=0;
307 while(nb_in_row<(A.mnel)[k2-1]){
308 dmumps_par->irn[k1]=k2; /* matrix indices start at 1 */
309 k1=k1+1;
310 nb_in_row=nb_in_row+1;
311 }
312 }
313 }
314
315 #if MUMPS_ARITH == MUMPS_ARITH_z
316 for(i=0;i<ne;i++){
317 ((dmumps_par->a)[i]).r = (A.R)[i];}
318 if(A.it == 1){
319 for(i=0;i<ne;i++){
320 ((dmumps_par->a)[i]).i = (A.I)[i];}
321 }else{
322 for(i=0;i<ne;i++){
323 ((dmumps_par->a)[i]).i = 0.0;}
324 }
325 #else
326 for(i=0;i<ne;i++){
327 ((dmumps_par->a)[i]) = (A.R)[i];}
328 #endif
329 dmumps_par->nz = ne;
330 }
331 else{
332 /* symmetric case */
333 k1=0;
334 i=0;
335 for (k2=1;k2<mA+1;k2++){
336 nb_in_row=0;
337 while(nb_in_row<(A.mnel)[k2-1]){
338 if( k2 >= (A.icol)[i]){
339 if(k1>=netrue){
340 Scierror(999,"%s: The matrix must be symmetric\n",function_name);
341 return 0;
342 }
343 (dmumps_par->jcn)[k1]=(A.icol)[i];
344 (dmumps_par->irn)[k1]=k2;
345 #if MUMPS_ARITH == MUMPS_ARITH_z
346 (dmumps_par->a)[k1].r=(A.R)[i];
347 if(A.it == 1){
348 ((dmumps_par->a)[k1]).i = (A.I)[i];}
349 else{
350 ((dmumps_par->a)[k1]).i = 0.0;}
351 #else
352 ((dmumps_par->a)[k1]) = (A.R)[i];
353 #endif
354 k1=k1+1;}
355
356 nb_in_row=nb_in_row+1;
357 i=i+1;
358 }
359 }
360 dmumps_par->nz = k1;
361 }
362
363 GetRhsVar(2,"i",&mjob,&njob,&ljob);
364 dmumps_par->job=*istk(ljob);
365
366 GetRhsVar(3,"i",&micntl,&nicntl,&licntl);
367 EXTRACT_FROM_SCILAB_TOARR(istk(licntl),dmumps_par->icntl,int,40);
368
369 GetRhsVar(4,"d",&mcntl,&ncntl,&lcntl);
370 EXTRACT_FROM_SCILAB_TOARR(stk(lcntl),dmumps_par->cntl,double,15);
371
372 GetRhsVar(5,"i",&mperm, &nperm, &lperm);
373 EXTRACT_FROM_SCILAB_TOPTR(IT_NOT_USED,istk(lperm),istk(lperm),(dmumps_par->perm_in),int,nA);
374
375 GetRhsCVar(6,"d",&it,&mcols,&ncols,&lcols,&licols);
376 EXTRACT_FROM_SCILAB_TOPTR(it,stk(lcols),stk(licols),(dmumps_par->colsca),double2,nA);
377
378 GetRhsCVar(7,"d",&it,&mrows,&nrows,&lrows,&lirows);
379 EXTRACT_FROM_SCILAB_TOPTR(it,stk(lrows),stk(lirows),(dmumps_par->rowsca),double2,nA);
380
381
382 /*
383 * To follow the "spirit" of the Matlab/Scilab interfaces, treat case of null
384 * space separately. In that case, we initialize lrhs and nrhs automatically,
385 * allocate the space needed, and do not rely on what is provided by the user
386 * in component RHS, that is not touched.
387 * Note that at the moment the user should not call the solution step combined
388 * with the factorization step when he/she sets icntl[25] to a non-zero value.
389 * Hence we suppose infog[28-1] is available and we can use it.
390 *
391 * For users of scilab/matlab, it would still be nice to be able to set ICNTL(25)=-1,
392 * and use JOB=6. If we want to make this functionality available, we should
393 * call separately job=2 and job=3 even if job=5 or 6 and set nrhs (and allocate
394 * space correctly) between job=2 and job=3 calls to MUMPS.
395 *
396 */
397 if ( dmumps_par->icntl[25-1] == -1 && dmumps_par->infog[28-1] > 0) {
398 dmumps_par->nrhs=dmumps_par->infog[28-1];
399 donullspace = dosolve;
400 }
401 else if ( dmumps_par->icntl[25-1] > 0 && dmumps_par->icntl[25-1] <= dmumps_par->infog[28-1] ) {
402 dmumps_par->nrhs=1;
403 donullspace = dosolve;
404 }
405 else {
406 donullspace=0;
407 }
408 if (donullspace) {
409 nRHS=dmumps_par->nrhs;
410 dmumps_par->lrhs=dmumps_par->n;
411 dmumps_par->rhs=(double2 *)malloc((dmumps_par->n)*(dmumps_par->nrhs)*sizeof(double2));
412 dmumps_par->icntl[19]=0;
413 }
414
415 else if(GetType(8)!=5){
416 /* Dense RHS */
417 GetRhsCVar(8,"d",&itRHS,&mRHS,&nRHS,&lRHS,&liRHS);
418
419 if((!dosolve) || (stk(lRHS)[0]) == -9999){
420 /* Could be dangerous ? See comment in Matlab interface */
421 EXTRACT_CMPLX_FROM_SCILAB_TOPTR(itRHS,stk(lRHS),stk(liRHS),(dmumps_par->rhs),double2,one);
422 }else{
423
424 dmumps_par->nrhs = nRHS;
425 dmumps_par->lrhs = mRHS;
426 if(mRHS!=nA){
427 Scierror(999,"%s: Incompatible number of rows in RHS\n",function_name);
428 }
429 dmumps_par->icntl[19]=0;
430 EXTRACT_CMPLX_FROM_SCILAB_TOPTR(itRHS,stk(lRHS),stk(liRHS),(dmumps_par->rhs),double2,(nRHS*mRHS));
431 }
432 }else{
433 /* Sparse RHS */
434 GetRhsVar(8,"s",&mRHS,&nRHS,&RHS_SPARSE);
435 dmumps_par->icntl[19]=1;
436 dmumps_par->nrhs = nRHS;
437 dmumps_par->lrhs = mRHS;
438 nz_RHS=RHS_SPARSE.nel;
439 dmumps_par->nz_rhs=nz_RHS;
440
441 irhs_ptr=(int*)malloc((nRHS+1)*sizeof(int));
442
443 dmumps_par->irhs_ptr=(int*)malloc((nRHS+1)*sizeof(int));
444 dmumps_par->irhs_sparse=(int*)malloc(nz_RHS*sizeof(int));
445 dmumps_par->rhs_sparse=(double2*)malloc(nz_RHS*sizeof(double2));
446 dmumps_par->rhs=(double2*)malloc((nRHS*mRHS)*sizeof(double2));
447 /* transform row-oriented sparse multiple rhs (scilab)
448 * into column-oriented sparse multiple rhs (mumps) */
449 k=0;
450 for(i=0;i<nRHS+1;i++){
451 irhs_ptr[i]=0;
452 dmumps_par->irhs_ptr[i]=0;}
453 for(i=1;i<mRHS+1;i++){
454 for(j=0;j<(RHS_SPARSE.mnel)[i-1];j++){
455 col_ind=(RHS_SPARSE.icol)[k];
456 k++;
457 ((dmumps_par->irhs_ptr)[col_ind])++;
458 }
459 }
460 (dmumps_par->irhs_ptr)[0]=1;
461 irhs_ptr[0]=(dmumps_par->irhs_ptr)[0];
462 for(i=1;i<nRHS+1;i++){
463 (dmumps_par->irhs_ptr)[i]=(dmumps_par->irhs_ptr)[i]+(dmumps_par->irhs_ptr)[i-1];
464 irhs_ptr[i]= (dmumps_par->irhs_ptr)[i];
465 }
466 k=RHS_SPARSE.nel-1;
467 for(i=mRHS;i>=1;i--){
468
469 for(j=0;j<(RHS_SPARSE.mnel)[i-1];j++){
470 col_ind=(RHS_SPARSE.icol)[k];
471 (dmumps_par->irhs_sparse)[irhs_ptr[col_ind]-2]=i;
472 #if MUMPS_ARITH == MUMPS_ARITH_z
473 if(RHS_SPARSE.it==1){
474 ((dmumps_par->rhs_sparse)[irhs_ptr[col_ind]-2]).r=RHS_SPARSE.R[k];
475 ((dmumps_par->rhs_sparse)[irhs_ptr[col_ind]-2]).i=RHS_SPARSE.I[k];
476 }else{
477 ((dmumps_par->rhs_sparse)[irhs_ptr[col_ind]-2]).r=RHS_SPARSE.R[k];
478 ((dmumps_par->rhs_sparse)[irhs_ptr[col_ind]-2]).i=0.0;
479 }
480 #else
481 (dmumps_par->rhs_sparse)[irhs_ptr[col_ind]-2]=RHS_SPARSE.R[k];
482 #endif
483 k--;
484 irhs_ptr[col_ind]=irhs_ptr[col_ind]-1;
485 }
486 }
487 MYFREE(irhs_ptr);
488 }
489
490 GetRhsVar(9,"i",&nv_schu,&mv_schu,&lv_schu);
491 dmumps_par-> size_schur=mv_schu;
492 EXTRACT_FROM_SCILAB_TOPTR(IT_NOT_USED,istk(lv_schu),istk(lv_schu),(dmumps_par->listvar_schur),int,dmumps_par->size_schur);
493 if(!dmumps_par->listvar_schur) dmumps_par->size_schur=0;
494
495 if(dmumps_par->size_schur > 0){
496 MYFREE(dmumps_par->schur);
497 if(!(dmumps_par->schur=(double2 *)malloc((dmumps_par->size_schur*dmumps_par->size_schur)*sizeof(double2)))){
498 Scierror(999,"%s: malloc Schur failed in intmumpsc.c\n",function_name);
499 }
500 dmumps_par->icntl[18]=1;
501 }else{
502 dmumps_par->icntl[18]=0;
503 }
504
505 /* Reduced RHS */
506 if ( dmumps_par->size_schur > 0 && dosolve ) {
507
508 if ( dmumps_par->icntl[26-1] == 2 ) {
509 /* REDRHS is on input */
510 GetRhsCVar(11,"d",&itREDRHS,&mREDRHS,&nREDRHS,&lREDRHS,&liREDRHS);
511 if (mREDRHS != dmumps_par->size_schur || nREDRHS != dmumps_par->nrhs ) {
512 Scierror(999,"%s: bad dimensions for REDRHS\n");
513 }
514 /* Fill dmumps_par->redrhs */
515 EXTRACT_CMPLX_FROM_SCILAB_TOPTR(itREDRHS,stk(lREDRHS),stk(liREDRHS),(dmumps_par->redrhs),double2,(nREDRHS*mREDRHS));
516 dmumps_par->lrhs=mREDRHS;
517 }
518
519 if ( dmumps_par->icntl[26-1] == 1 ) {
520 /* REDRHS on output. Must be allocated before the call */
521 MYFREE(dmumps_par->redrhs);
522 if(!(dmumps_par->redrhs=(double2 *)malloc((dmumps_par->size_schur*dmumps_par->nrhs)*sizeof(double2)))){
523 Scierror(999,"%s: malloc redrhs failed in intmumpsc.c\n",function_name);
524 }
525 }
526 }
527
528 /* call C interface to MUMPS */
529 dmumps_c(dmumps_par);
530
531 }
532 }
533
534 if(*istk(ljob)==-2){
535 return 0;
536 }else{
537
538 CheckLhs(11,11);
539
540 EXTRACT_INT_FROM_C_TO_SCILAB(1,linfog,(dmumps_par->infog),one,temp1,one);
541
542 EXTRACT_DOUBLE_FROM_C_TO_SCILAB(2,it,lrinfog,lrinfog,(dmumps_par->rinfog),one,temp2,one);
543
544 if(dmumps_par->rhs && dosolve){ /* Just to know if solution step was called */
545 it =1;
546 EXTRACT_CMPLX_FROM_C_TO_SCILAB(3,it,lrhsout,lrhsouti,(dmumps_par->rhs),nA,nRHS,one);
547
548 }else{
549 it=1;
550 EXTRACT_CMPLX_FROM_C_TO_SCILAB(3,it,lrhsout,lrhsouti,(dmumps_par->rhs),one,one,one);
551 }
552
553 ptr_int = (int *)dmumps_par;
554 inst_address = (int) ptr_int;
555 EXTRACT_INT_FROM_C_TO_SCILAB(4,linstout,&inst_address,one,one,one);
556
557
558 temp4=dmumps_par->size_schur;
559 if(temp4>0){
560 it=1;
561 EXTRACT_CMPLX_FROM_C_TO_SCILAB(5,it,lschurout,lschurouti,(dmumps_par->schur),temp4,temp4,one);
562 }else{
563 it=1;
564 EXTRACT_CMPLX_FROM_C_TO_SCILAB(5,it,lschurout,lschurouti,(dmumps_par->schur),one,one,one);
565 }
566
567 /* REDRHS on output */
568 it=1;
569 if ( dmumps_par->icntl[26-1]==1 && dmumps_par->size_schur > 0 && dosolve ) {
570 len1=dmumps_par->size_schur;
571 len2=dmumps_par->nrhs;
572 }
573 else {
574 len1=1;
575 len2=1;
576 }
577 it=1;
578 EXTRACT_CMPLX_FROM_C_TO_SCILAB(6,it,stkptr,stkptri,(dmumps_par->redrhs),len1,len2,one)
579
580
581 MYFREE(dmumps_par->redrhs);
582 MYFREE(dmumps_par->schur);
583 MYFREE(dmumps_par->irhs_ptr);
584 MYFREE(dmumps_par->irhs_sparse);
585 MYFREE(dmumps_par->rhs_sparse);
586 MYFREE(dmumps_par->rhs);
587
588 /* temp3=dmumps_par->deficiency;*/
589 temp3=dmumps_par->infog[27];
590 EXTRACT_INT_FROM_C_TO_SCILAB(7,lpivnul_list,(dmumps_par->pivnul_list),one,temp3,one);
591
592 EXTRACT_INT_FROM_C_TO_SCILAB(8,lsymperm,(dmumps_par->sym_perm),one,nA,one);
593
594 EXTRACT_INT_FROM_C_TO_SCILAB(9,lunsperm,(dmumps_par->uns_perm),one,nA,one);
595
596 nicntl=40;
597 EXTRACT_INT_FROM_C_TO_SCILAB(10,licntl,(dmumps_par->icntl),one,nicntl,one);
598 ncntl=15;
599 EXTRACT_DOUBLE_FROM_C_TO_SCILAB(11,it,lcntl,lcntl,(dmumps_par->cntl),one,ncntl,one);
600 return 0;
601
602 }
603 }
604
605
606 static GenericTable Tab[]={
607 #if MUMPS_ARITH == MUMPS_ARITH_z
608 {(Myinterfun) sci_gateway, dmumpsc,"zmumpsc"}
609 #else
610 {(Myinterfun) sci_gateway, dmumpsc,"dmumpsc"}
611 #endif
612 };
613
614 #if MUMPS_ARITH == MUMPS_ARITH_z
C2F(scizmumps)615 int C2F(scizmumps)()
616 #else
617 int C2F(scidmumps)()
618 #endif
619 {Rhs = Max(0, Rhs);
620 (*(Tab[Fin-1].f))(Tab[Fin-1].name,Tab[Fin-1].F);
621 return 0;
622 }
623