181 {
182
183
184
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
193 int linfog, lrinfog, lrhsout,lrhsouti, linstout, lschurout, lschurouti, ldef;
194 int lpivnul_list, lmapp, lsymperm, lunsperm;
195 int one=1, temp1=80, temp2=40, temp3, temp4;
196 int it, itRHS, itREDRHS;
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
207 int len1, len2;
208
209 int stkptr, stkptri;
210
211
212 int inst_address;
213 int ne,inst;
214 int *irn_in,*jcn_in;
215
216
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;
231
232 int dosolve=0;
233 int donullspace=0;
234 int doanal = 0;
235
236 CheckRhs(11,12);
237
238
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
245 GetRhsVar(1,"i",&msym,&nsym,&lsym);
246 dmumps_par->
sym=*istk(lsym);
247 dmumps_par->
job = -1;
252 it=1;
253 }else{
254
255 GetRhsVar(10,"i",&mint,&nint,&lint);
256 inst_address=*istk(lint);
257 ptr_int = (int *) inst_address;
258
260 if(*istk(ljob) == -2){
261 dmumps_par->
job = -2;
264 }else{
265
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;
275
276 if(dmumps_par->
sym != 0){
277 netrue = (nA+ne)/2;
278 }else{
279 netrue = ne;
280 }
281
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));
291 }
292
293
294 if ((dmumps_par->
sym)==0){
295
296
297
298
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;
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
330 }
331 else{
332
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 }
361 }
362
363 GetRhsVar(2,"i",&mjob,&njob,&ljob);
364 dmumps_par->
job=*istk(ljob);
365
366 GetRhsVar(3,"i",&micntl,&nicntl,&licntl);
368
369 GetRhsVar(4,"d",&mcntl,&ncntl,&lcntl);
371
372 GetRhsVar(5,"i",&mperm, &nperm, &lperm);
374
375 GetRhsCVar(6,"d",&it,&mcols,&ncols,&lcols,&licols);
377
378 GetRhsCVar(7,"d",&it,&mrows,&nrows,&lrows,&lirows);
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
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] ) {
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
417 GetRhsCVar(8,"d",&itRHS,&mRHS,&nRHS,&lRHS,&liRHS);
418
419 if((!dosolve) || (stk(lRHS)[0]) == -9999){
420
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
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
448
449 k=0;
450 for(i=0;i<nRHS+1;i++){
451 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 }
461 irhs_ptr[0]=(dmumps_par->
irhs_ptr)[0];
462 for(i=1;i<nRHS+1;i++){
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];
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 }
488 }
489
490 GetRhsVar(9,"i",&nv_schu,&mv_schu,&lv_schu);
491 dmumps_par-> size_schur=mv_schu;
494
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
506 if ( dmumps_par->
size_schur > 0 && dosolve ) {
507
508 if ( dmumps_par->
icntl[26-1] == 2 ) {
509
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
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
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
530
531 }
532 }
533
534 if(*istk(ljob)==-2){
535 return 0;
536 }else{
537
538 CheckLhs(11,11);
539
541
543
544 if(dmumps_par->
rhs && dosolve){
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;
556
557
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
568 it=1;
569 if ( dmumps_par->
icntl[26-1]==1 && dmumps_par->
size_schur > 0 && dosolve ) {
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
587
588
589 temp3=dmumps_par->
infog[27];
591
593
595
596 nicntl=60;
598 ncntl=15;
600 return 0;
601
602 }
603}
void MUMPS_CALL dmumps_c(DMUMPS_STRUC_C *dmumps_par)
#define EXTRACT_FROM_SCILAB_TOARR(ptr_scilab, mumpsarray, type, length)
void DMUMPS_free(DMUMPS_STRUC_C **dmumps_par)
#define EXTRACT_DOUBLE_FROM_C_TO_SCILAB(num, it, ptr_scilab1, ptr_scilab2, mumpsptr, length1, length2, one)
#define EXTRACT_FROM_SCILAB_TOPTR(it, ptr_scilab1, ptr_scilab2, mumpspointer, type, length)
void DMUMPS_alloc(DMUMPS_STRUC_C **dmumps_par)
#define EXTRACT_INT_FROM_C_TO_SCILAB(num, ptr_scilab, mumpsptr, length1, length2, one)
DMUMPS_COMPLEX * rhs_sparse
MUMPS_INT * listvar_schur