OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mumpsmex.c File Reference
#include "mex.h"

Go to the source code of this file.

Macros

#define MUMPS_ARITH_d   2
#define MUMPS_ARITH_z   8
#define SYM   (prhs[0])
#define JOB   (prhs[1])
#define ICNTL_IN   (prhs[2])
#define CNTL_IN   (prhs[3])
#define PERM_IN   (prhs[4])
#define COLSCA_IN   (prhs[5])
#define ROWSCA_IN   (prhs[6])
#define RHS_IN   (prhs[7])
#define VAR_SCHUR   (prhs[8])
#define INST   (prhs[9])
#define REDRHS_IN   (prhs[10])
#define KEEP_IN   (prhs[11])
#define DKEEP_IN   (prhs[12])
#define A_IN   (prhs[13])
#define INFO_OUT   (plhs[0])
#define RINFO_OUT   (plhs[1])
#define RHS_OUT   (plhs[2])
#define INST_OUT   (plhs[3])
#define SCHUR_OUT   (plhs[4])
#define REDRHS_OUT   (plhs[5])
#define PIVNUL_LIST   (plhs[6])
#define PERM_OUT   (plhs[7])
#define UNS_PERM   (plhs[8])
#define ICNTL_OUT   (plhs[9])
#define CNTL_OUT   (plhs[10])
#define COLSCA_OUT   (plhs[11])
#define ROWSCA_OUT   (plhs[12])
#define KEEP_OUT   (plhs[13])
#define DKEEP_OUT   (plhs[14])
#define MYMALLOC(ptr, l, type)
#define MYFREE(ptr)
#define EXTRACT_FROM_MATLAB_TOPTR(mxcomponent, mumpspointer, type, length)
#define EXTRACT_SCALING_FROM_MATLAB_TOPTR(mxcomponent, mumpspointer, is_a_pointer_from_mumps, length)
#define EXTRACT_FROM_MATLAB_TOARR(mxcomponent, mumpsarray, type, length)
#define EXTRACT_FROM_MATLAB_TOVAL(mxcomponent, mumpsvalue)
#define EXTRACT_FROM_C_TO_MATLAB(mxcomponent, mumpspointer, length)

Functions

void DMUMPS_free (DMUMPS_STRUC_C **dmumps_par)
void DMUMPS_alloc (DMUMPS_STRUC_C **dmumps_par)
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

Macro Definition Documentation

◆ A_IN

#define A_IN   (prhs[13])

Definition at line 45 of file mumpsmex.c.

◆ CNTL_IN

#define CNTL_IN   (prhs[3])

Definition at line 35 of file mumpsmex.c.

◆ CNTL_OUT

#define CNTL_OUT   (plhs[10])

Definition at line 57 of file mumpsmex.c.

◆ COLSCA_IN

#define COLSCA_IN   (prhs[5])

Definition at line 37 of file mumpsmex.c.

◆ COLSCA_OUT

#define COLSCA_OUT   (plhs[11])

Definition at line 58 of file mumpsmex.c.

◆ DKEEP_IN

#define DKEEP_IN   (prhs[12])

Definition at line 44 of file mumpsmex.c.

◆ DKEEP_OUT

#define DKEEP_OUT   (plhs[14])

Definition at line 61 of file mumpsmex.c.

◆ EXTRACT_FROM_C_TO_MATLAB

#define EXTRACT_FROM_C_TO_MATLAB ( mxcomponent,
mumpspointer,
length )
Value:
if(mumpspointer == 0){ \
mxcomponent = mxCreateDoubleMatrix (1, 1, mxREAL); \
ptr_matlab = mxGetPr (mxcomponent); \
ptr_matlab[0] = -9999; \
}else{ \
mxcomponent = mxCreateDoubleMatrix (1,length,mxREAL); \
ptr_matlab = mxGetPr (mxcomponent); \
for(i=0;i<length;i++){ \
ptr_matlab[i]=(double)(mumpspointer)[i]; \
} \
}

Definition at line 117 of file mumpsmex.c.

117#define EXTRACT_FROM_C_TO_MATLAB(mxcomponent,mumpspointer,length) \
118 if(mumpspointer == 0){ \
119 mxcomponent = mxCreateDoubleMatrix (1, 1, mxREAL); \
120 ptr_matlab = mxGetPr (mxcomponent); \
121 ptr_matlab[0] = -9999; \
122 }else{ \
123 mxcomponent = mxCreateDoubleMatrix (1,length,mxREAL); \
124 ptr_matlab = mxGetPr (mxcomponent); \
125 for(i=0;i<length;i++){ \
126 ptr_matlab[i]=(double)(mumpspointer)[i]; \
127 } \
128 }

◆ EXTRACT_FROM_MATLAB_TOARR

#define EXTRACT_FROM_MATLAB_TOARR ( mxcomponent,
mumpsarray,
type,
length )
Value:
ptr_matlab = mxGetPr(mxcomponent); \
if(ptr_matlab[0] != -9999){ \
for(i=0;i<length;i++){ \
if(ptr_matlab[i] != -9998){ \
mumpsarray[i] = ptr_matlab[i]; \
} \
} \
}

Definition at line 101 of file mumpsmex.c.

101#define EXTRACT_FROM_MATLAB_TOARR(mxcomponent,mumpsarray,type,length) \
102 ptr_matlab = mxGetPr(mxcomponent); \
103 if(ptr_matlab[0] != -9999){ \
104 for(i=0;i<length;i++){ \
105 if(ptr_matlab[i] != -9998){ \
106 mumpsarray[i] = ptr_matlab[i]; \
107 } \
108 } \
109 }

◆ EXTRACT_FROM_MATLAB_TOPTR

#define EXTRACT_FROM_MATLAB_TOPTR ( mxcomponent,
mumpspointer,
type,
length )
Value:
ptr_matlab = mxGetPr(mxcomponent); \
if(ptr_matlab[0] != -9999){ \
MYFREE(mumpspointer); \
MYMALLOC(mumpspointer,length,type); \
for(i=0;i<length;i++){ \
mumpspointer[i] = ptr_matlab[i]; \
} \
}

Definition at line 74 of file mumpsmex.c.

74#define EXTRACT_FROM_MATLAB_TOPTR(mxcomponent,mumpspointer,type,length) \
75 ptr_matlab = mxGetPr(mxcomponent); \
76 if(ptr_matlab[0] != -9999){ \
77 MYFREE(mumpspointer); \
78 MYMALLOC(mumpspointer,length,type); \
79 for(i=0;i<length;i++){ \
80 mumpspointer[i] = ptr_matlab[i]; \
81 } \
82 }

◆ EXTRACT_FROM_MATLAB_TOVAL

#define EXTRACT_FROM_MATLAB_TOVAL ( mxcomponent,
mumpsvalue )
Value:
ptr_matlab = mxGetPr(mxcomponent); \
if(ptr_matlab[0] != -9999){ \
mumpsvalue = ptr_matlab[0]; \
}

Definition at line 111 of file mumpsmex.c.

111#define EXTRACT_FROM_MATLAB_TOVAL(mxcomponent,mumpsvalue) \
112 ptr_matlab = mxGetPr(mxcomponent); \
113 if(ptr_matlab[0] != -9999){ \
114 mumpsvalue = ptr_matlab[0]; \
115 }

◆ EXTRACT_SCALING_FROM_MATLAB_TOPTR

#define EXTRACT_SCALING_FROM_MATLAB_TOPTR ( mxcomponent,
mumpspointer,
is_a_pointer_from_mumps,
length )
Value:
ptr_matlab = mxGetPr(mxcomponent); \
if( ptr_matlab[0] != -9999 && ! (is_a_pointer_from_mumps) ) { \
MYFREE(mumpspointer); \
MYMALLOC(mumpspointer,length,double); \
for(i=0;i<length;i++){ \
mumpspointer[i] = ptr_matlab[i]; \
} \
}

Definition at line 91 of file mumpsmex.c.

91#define EXTRACT_SCALING_FROM_MATLAB_TOPTR(mxcomponent,mumpspointer,is_a_pointer_from_mumps,length) \
92 ptr_matlab = mxGetPr(mxcomponent); \
93 if( ptr_matlab[0] != -9999 && ! (is_a_pointer_from_mumps) ) { \
94 MYFREE(mumpspointer); \
95 MYMALLOC(mumpspointer,length,double); \
96 for(i=0;i<length;i++){ \
97 mumpspointer[i] = ptr_matlab[i]; \
98 } \
99 }

◆ ICNTL_IN

#define ICNTL_IN   (prhs[2])

Definition at line 34 of file mumpsmex.c.

◆ ICNTL_OUT

#define ICNTL_OUT   (plhs[9])

Definition at line 56 of file mumpsmex.c.

◆ INFO_OUT

#define INFO_OUT   (plhs[0])

Definition at line 47 of file mumpsmex.c.

◆ INST

id INST   (prhs[9])

Definition at line 41 of file mumpsmex.c.

◆ INST_OUT

#define INST_OUT   (plhs[3])

Definition at line 50 of file mumpsmex.c.

◆ JOB

id JOB   (prhs[1])

Definition at line 33 of file mumpsmex.c.

◆ KEEP_IN

#define KEEP_IN   (prhs[11])

Definition at line 43 of file mumpsmex.c.

◆ KEEP_OUT

#define KEEP_OUT   (plhs[13])

Definition at line 60 of file mumpsmex.c.

◆ MUMPS_ARITH_d

#define MUMPS_ARITH_d   2

Definition at line 3 of file mumpsmex.c.

◆ MUMPS_ARITH_z

#define MUMPS_ARITH_z   8

Definition at line 4 of file mumpsmex.c.

◆ MYFREE

#define MYFREE ( ptr)
Value:
if(ptr){ \
free(ptr); \
ptr = 0; \
}

Definition at line 68 of file mumpsmex.c.

68#define MYFREE(ptr) \
69 if(ptr){ \
70 free(ptr); \
71 ptr = 0; \
72 }

◆ MYMALLOC

#define MYMALLOC ( ptr,
l,
type )
Value:
if(!(ptr = (type *) malloc(l*sizeof(type)))){ \
mexErrMsgTxt ("Malloc failed in mumpsmex.c"); \
}

Definition at line 63 of file mumpsmex.c.

63#define MYMALLOC(ptr,l,type) \
64 if(!(ptr = (type *) malloc(l*sizeof(type)))){ \
65 mexErrMsgTxt ("Malloc failed in mumpsmex.c"); \
66 }

◆ PERM_IN

#define PERM_IN   (prhs[4])

Definition at line 36 of file mumpsmex.c.

◆ PERM_OUT

#define PERM_OUT   (plhs[7])

Definition at line 54 of file mumpsmex.c.

◆ PIVNUL_LIST

id PIVNUL_LIST   (plhs[6])

Definition at line 53 of file mumpsmex.c.

◆ REDRHS_IN

#define REDRHS_IN   (prhs[10])

Definition at line 42 of file mumpsmex.c.

◆ REDRHS_OUT

#define REDRHS_OUT   (plhs[5])

Definition at line 52 of file mumpsmex.c.

◆ RHS_IN

#define RHS_IN   (prhs[7])

Definition at line 39 of file mumpsmex.c.

◆ RHS_OUT

#define RHS_OUT   (plhs[2])

Definition at line 49 of file mumpsmex.c.

◆ RINFO_OUT

#define RINFO_OUT   (plhs[1])

Definition at line 48 of file mumpsmex.c.

◆ ROWSCA_IN

#define ROWSCA_IN   (prhs[6])

Definition at line 38 of file mumpsmex.c.

◆ ROWSCA_OUT

#define ROWSCA_OUT   (plhs[12])

Definition at line 59 of file mumpsmex.c.

◆ SCHUR_OUT

#define SCHUR_OUT   (plhs[4])

Definition at line 51 of file mumpsmex.c.

◆ SYM

#define SYM   (prhs[0])

Definition at line 32 of file mumpsmex.c.

◆ UNS_PERM

id UNS_PERM   (plhs[8])

Definition at line 55 of file mumpsmex.c.

◆ VAR_SCHUR

#define VAR_SCHUR   (prhs[8])

Definition at line 40 of file mumpsmex.c.

Function Documentation

◆ DMUMPS_alloc()

void DMUMPS_alloc ( DMUMPS_STRUC_C ** dmumps_par)

Definition at line 204 of file mumpsmex.c.

204 {
205
206 MYMALLOC((*dmumps_par),1,DMUMPS_STRUC_C);
207 (*dmumps_par)->irn = NULL;
208 (*dmumps_par)->jcn = NULL;
209 (*dmumps_par)->a = NULL;
210 (*dmumps_par)->irn_loc = NULL;
211 (*dmumps_par)->jcn_loc = NULL;
212 (*dmumps_par)->a_loc = NULL;
213 (*dmumps_par)->eltptr = NULL;
214 (*dmumps_par)->eltvar = NULL;
215 (*dmumps_par)->a_elt = NULL;
216 (*dmumps_par)->perm_in = NULL;
217 (*dmumps_par)->colsca = NULL;
218 (*dmumps_par)->rowsca = NULL;
219 (*dmumps_par)->rhs = NULL;
220 (*dmumps_par)->redrhs = NULL;
221 (*dmumps_par)->rhs_sparse = NULL;
222 (*dmumps_par)->irhs_sparse = NULL;
223 (*dmumps_par)->irhs_ptr = NULL;
224 (*dmumps_par)->pivnul_list = NULL;
225 (*dmumps_par)->listvar_schur = NULL;
226 (*dmumps_par)->schur = NULL;
227 (*dmumps_par)->sym_perm = NULL;
228 (*dmumps_par)->uns_perm = NULL;
229}
MYMALLOC
#define MYMALLOC(ptr, l, type)
Definition mumpsmex.c:63
DMUMPS_STRUC_C
Definition dmumps_c.h:42

◆ DMUMPS_free()

void DMUMPS_free ( DMUMPS_STRUC_C ** dmumps_par)

Definition at line 174 of file mumpsmex.c.

174 {
175 if(*dmumps_par){
176 MYFREE( (*dmumps_par)->irn );
177 MYFREE( (*dmumps_par)->jcn );
178 MYFREE( (*dmumps_par)->a );
179 MYFREE( (*dmumps_par)->irn_loc );
180 MYFREE( (*dmumps_par)->jcn_loc );
181 MYFREE( (*dmumps_par)->a_loc );
182 MYFREE( (*dmumps_par)->eltptr );
183 MYFREE( (*dmumps_par)->eltvar );
184 MYFREE( (*dmumps_par)->a_elt );
185 MYFREE( (*dmumps_par)->perm_in );
186 /* colsca/rowsca might have been allocated by
187 * MUMPS but in that case the corresponding pointer
188 * is already equal to 0 so line below will do nothing */
189 MYFREE( (*dmumps_par)->colsca );
190 MYFREE( (*dmumps_par)->rowsca );
191 MYFREE( (*dmumps_par)->pivnul_list );
192 MYFREE( (*dmumps_par)->listvar_schur );
193 MYFREE( (*dmumps_par)->sym_perm );
194 MYFREE( (*dmumps_par)->uns_perm );
195 MYFREE( (*dmumps_par)->irhs_ptr);
196 MYFREE( (*dmumps_par)->irhs_sparse);
197 MYFREE( (*dmumps_par)->rhs_sparse);
198 MYFREE( (*dmumps_par)->rhs);
199 MYFREE( (*dmumps_par)->redrhs);
200 MYFREE(*dmumps_par);
201 }
202}
MYFREE
#define MYFREE(ptr)
Definition mumpsmex.c:68

◆ mexFunction()

void mexFunction ( int nlhs,
mxArray * plhs[],
int nrhs,
const mxArray * prhs[] )

Definition at line 231 of file mumpsmex.c.

232 {
233
234 int i,j,pos;
235 int *ptr_int;
236 double *ptr_matlab;
237#if MUMPS_ARITH == MUMPS_ARITH_z
238 double *ptri_matlab;
239#endif
240 mwSize tmp_m,tmp_n;
241
242 /* C pointer for input parameters */
243 size_t inst_address;
244 mwSize n,m,ne, netrue ;
245 int job;
246 mwIndex *irn_in,*jcn_in;
247
248 /* variable for multiple and sparse rhs */
249 int posrhs;
250 mwSize nbrhs,ldrhs, nz_rhs;
251 mwIndex *irhs_ptr, *irhs_sparse;
252 double *rhs_sparse;
253#if MUMPS_ARITH == MUMPS_ARITH_z
254 double *im_rhs_sparse;
255#endif
256
257 DMUMPS_STRUC_C *dmumps_par;
258 int dosolve = 0;
259 int donullspace = 0;
260 int doanalysis = 0;
261 int dofactorize = 0;
262
263
264 EXTRACT_FROM_MATLAB_TOVAL(JOB,job);
265
266 doanalysis = (job == 1 || job == 4 || job == 6);
267 dofactorize = (job == 2 || job == 4 || job == 5 || job == 6);
268 dosolve = (job == 3 || job == 5 || job == 6);
269
270 if(job == -1){
271 DMUMPS_alloc(&dmumps_par);
272 EXTRACT_FROM_MATLAB_TOVAL(SYM,dmumps_par->sym);
273 dmumps_par->job = -1;
274 dmumps_par->par = 1;
275 dmumps_c(dmumps_par);
276 dmumps_par->nz = -1;
277 dmumps_par->nz_alloc = -1;
278 }else{
279 EXTRACT_FROM_MATLAB_TOVAL(INST,inst_address);
280 ptr_int = (int *) inst_address;
281
282 dmumps_par = (DMUMPS_STRUC_C *) ptr_int;
283
284 if(job == -2){
285 dmumps_par->job = -2;
286 dmumps_c(dmumps_par);
287 /* If colsca/rowsca were freed by MUMPS,
288 dmumps_par->colsca/rowsca are now null.
289 Application of MYFREE in call below thus ok */
290 DMUMPS_free(&dmumps_par);
291 }else{
292
293 /* check of input arguments */
294 n = mxGetN(A_IN);
295 m = mxGetM(A_IN);
296
297 if (!mxIsSparse(A_IN) || n != m )
298 mexErrMsgTxt("Input matrix must be a sparse square matrix");
299
300 jcn_in = mxGetJc(A_IN);
301 ne = jcn_in[n];
302 irn_in = mxGetIr(A_IN);
303 dmumps_par->n = (int)n;
304 if(dmumps_par->n != n)
305 mexErrMsgTxt("Input is too big; will not work...barfing out\n");
306
307 if(dmumps_par->sym != 0)
308 netrue = (n+ne)/2;
309 else
310 netrue = ne;
311
312 if(dmumps_par->nz_alloc < netrue || dmumps_par->nz_alloc >= 2*netrue){
313 MYFREE(dmumps_par->jcn);
314 MYFREE(dmumps_par->irn);
315 MYFREE(dmumps_par->a);
316 MYMALLOC((dmumps_par->jcn),(int)netrue,int);
317 MYMALLOC((dmumps_par->irn),(int)netrue,int);
318 MYMALLOC((dmumps_par->a),(int)netrue,double2);
319 dmumps_par->nz_alloc = (int)netrue;
320 if (dmumps_par->nz_alloc != netrue)
321 mexErrMsgTxt("Input is too big; will not work...barfing out\n");
322 }
323
324
325 if(dmumps_par->sym == 0){
326 /* if analysis already performed then we only need to read
327 numerical values
328 Note that we suppose that matlab did not change the internal
329 format of the matrix between the 2 calls */
330 if(doanalysis){
331 /* || dmumps_par->info[22] == 0 */
332 for(i=0;i<dmumps_par->n;i++){
333 for(j=jcn_in[i];j<jcn_in[i+1];j++){
334 (dmumps_par->jcn)[j] = i+1;
335 (dmumps_par->irn)[j] = ((int)irn_in[j])+1;
336 }
337 }
338 }
339 dmumps_par->nz = (int)ne;
340 if( dmumps_par->nz != ne)
341 mexErrMsgTxt("Input is too big; will not work...barfing out\n");
342#if MUMPS_ARITH == MUMPS_ARITH_z
343 ptr_matlab = mxGetPr(A_IN);
344 for(i=0;i<dmumps_par->nz;i++){
345 ((dmumps_par->a)[i]).r = ptr_matlab[i];
346 }
347 ptr_matlab = mxGetPi(A_IN);
348 if(ptr_matlab){
349 for(i=0;i<dmumps_par->nz;i++){
350 ((dmumps_par->a)[i]).i = ptr_matlab[i];
351 }
352 }else{
353 for(i=0;i<dmumps_par->nz;i++){
354 ((dmumps_par->a)[i]).i = 0.0;
355 }
356 }
357#else
358 ptr_matlab = mxGetPr(A_IN);
359 for(i=0;i<dmumps_par->nz;i++){
360 (dmumps_par->a)[i] = ptr_matlab[i];
361 }
362#endif
363 }else{
364 /* in the symmetric case we do not need to check doanalysis */
365 pos = 0;
366 ptr_matlab = mxGetPr(A_IN);
367#if MUMPS_ARITH == MUMPS_ARITH_z
368 ptri_matlab = mxGetPi(A_IN);
369#endif
370 for(i=0;i<dmumps_par->n;i++){
371 for(j=jcn_in[i];j<jcn_in[i+1];j++){
372 if(irn_in[j] >= i){
373 if(pos >= netrue)
374 mexErrMsgTxt("Input matrix must be symmetric");
375 (dmumps_par->jcn)[pos] = i+1;
376 (dmumps_par->irn)[pos] = (int)irn_in[j]+1;
377#if MUMPS_ARITH == MUMPS_ARITH_z
378 ((dmumps_par->a)[pos]).r = ptr_matlab[j];
379 if(ptri_matlab){
380 ((dmumps_par->a)[pos]).i = ptri_matlab[j];
381 }else{
382 ((dmumps_par->a)[pos]).i = 0.0;
383 }
384#else
385 (dmumps_par->a)[pos] = ptr_matlab[j];
386#endif
387 pos++;
388 }
389 }
390 }
391 dmumps_par->nz = pos;
392 }
393
394
395 EXTRACT_FROM_MATLAB_TOVAL(JOB,dmumps_par->job);
396 EXTRACT_FROM_MATLAB_TOARR(ICNTL_IN,dmumps_par->icntl,int,60);
397 EXTRACT_FROM_MATLAB_TOARR(CNTL_IN,dmumps_par->cntl,double,15);
398 EXTRACT_FROM_MATLAB_TOPTR(PERM_IN,(dmumps_par->perm_in),int,((int)n));
399
400 /* colsca and rowsca are treated differently: it may happen that
401 dmumps_par-> colsca is nonzero because it was set to a nonzero
402 value on output (COLSCA_OUT) from MUMPS. Unfortunately if scaling
403 was on output, one cannot currently provide scaling on input
404 afterwards without reinitializing the instance */
405
406 EXTRACT_SCALING_FROM_MATLAB_TOPTR(COLSCA_IN,(dmumps_par->colsca),(dmumps_par->colsca_from_mumps),((int)n)); /* type always double */
407 EXTRACT_SCALING_FROM_MATLAB_TOPTR(ROWSCA_IN,(dmumps_par->rowsca),(dmumps_par->rowsca_from_mumps),((int)n)); /* type always double */
408
409 EXTRACT_FROM_MATLAB_TOARR(KEEP_IN,dmumps_par->keep,int,500);
410 EXTRACT_FROM_MATLAB_TOARR(DKEEP_IN,dmumps_par->dkeep,double,230);
411
412 dmumps_par->size_schur = (int)mxGetN(VAR_SCHUR);
413 EXTRACT_FROM_MATLAB_TOPTR(VAR_SCHUR,(dmumps_par->listvar_schur),int,dmumps_par->size_schur);
414 if(!dmumps_par->listvar_schur) dmumps_par->size_schur = 0;
415
416 ptr_matlab = mxGetPr (RHS_IN);
417
418/*
419 * To follow the "spirit" of the Matlab/Scilab interfaces, treat case of null
420 * space separately. In that case, we initialize lrhs and nrhs, automatically
421 * allocate the space needed, and do not rely on what is provided by the user
422 * in component RHS, that is not touched.
423 *
424 * Note that, at the moment, the user should not call the solution step combined
425 * with the factorization step when he/she sets icntl[25-1] to a non-zero value.
426 * Hence we suppose in the following that infog[28-1] is available and that we
427 * can use it.
428 *
429 * For users of scilab/matlab, it would still be nice to be able to set ICNTL(25)=-1,
430 * and use JOB=6. If we want to make such a feature available, we should
431 * call separately job=2 and job=3 even if job=5 or 6 and set nbrhs (and allocate
432 * space correctly) between job=2 and job=3 calls to MUMPS.
433 *
434 */
435 if ( dmumps_par->icntl[25-1] == -1 && dmumps_par->infog[28-1] > 0 ) {
436 dmumps_par->nrhs=dmumps_par->infog[28-1];
437 donullspace = dosolve;
438 }
439 else if ( dmumps_par->icntl[25-1] > 0 && dmumps_par->icntl[25-1] <= dmumps_par->infog[28-1] ) {
440 dmumps_par->nrhs=1;
441 donullspace = dosolve;
442 }
443 else {
444 donullspace=0;
445 }
446 if (donullspace) {
447 nbrhs=dmumps_par->nrhs; ldrhs=n;
448 dmumps_par->lrhs=(int)n;
449 MYMALLOC((dmumps_par->rhs),((dmumps_par->n)*(dmumps_par->nrhs)),double2);
450 }
451 else if((!dosolve) || ptr_matlab[0] == -9999 ) { /* rhs not already provided, or not used */
452/* Case where dosolve is true and ptr_matlab[0]=-9999, this could cause problems:
453 * 1/ RHS was not initialized while it should have been
454 * 2/ RHS was explicitely initialized to -9999 but is not allocated of the right size
455 */
456 EXTRACT_CMPLX_FROM_MATLAB_TOPTR(RHS_IN,(dmumps_par->rhs),double,1);
457 }else{
458 nbrhs = mxGetN(RHS_IN);
459 ldrhs = mxGetM(RHS_IN);
460 dmumps_par->nrhs = (int)nbrhs;
461 dmumps_par->lrhs = (int)ldrhs;
462 if(ldrhs != n){
463 mexErrMsgTxt ("Incompatible number of rows in RHS");
464 }
465 if (!mxIsSparse(RHS_IN)){ /* full rhs */
466 dmumps_par->icntl[20-1] = 0;
467 EXTRACT_CMPLX_FROM_MATLAB_TOPTR(RHS_IN,(dmumps_par->rhs),double,(int)( dmumps_par->nrhs*ldrhs));
468 }else{ /* sparse rhs */
469 /* printf("sparse RHS ldrhs = %d nrhs = %d\n",ldrhs,nbrhs); */
470 if (dmumps_par->icntl[30-1] == 0) {
471 /* A-1 feature was not requested => we are in the standard
472 * sparse RHS case and thus we set ICNTL(20) accordingly. */
473 dmumps_par->icntl[20-1] = 1;
474 }
475 irhs_ptr = mxGetJc(RHS_IN);
476 irhs_sparse = mxGetIr(RHS_IN);
477 rhs_sparse = mxGetPr(RHS_IN);
478#if MUMPS_ARITH == MUMPS_ARITH_z
479 im_rhs_sparse = mxGetPi(RHS_IN);
480#endif
481 nz_rhs = irhs_ptr[nbrhs];
482 dmumps_par->nz_rhs = (int)nz_rhs;
483
484 MYMALLOC((dmumps_par->irhs_ptr),(dmumps_par->nrhs+1),int);
485 MYMALLOC((dmumps_par->irhs_sparse), dmumps_par->nz_rhs,int);
486 MYMALLOC((dmumps_par->rhs_sparse), dmumps_par->nz_rhs,double2);
487 /* dmumps_par->rhs will store the solution*/
488 MYMALLOC((dmumps_par->rhs),((dmumps_par->nrhs*dmumps_par->lrhs)),double2);
489
490 for(i=0;i< dmumps_par->nrhs;i++){
491 for(j=irhs_ptr[i];j<irhs_ptr[i+1];j++){
492 (dmumps_par->irhs_sparse)[j] = irhs_sparse[j]+1;
493 }
494 (dmumps_par->irhs_ptr)[i] = irhs_ptr[i]+1;
495 }
496 (dmumps_par->irhs_ptr)[dmumps_par->nrhs] = dmumps_par->nz_rhs+1;
497#if MUMPS_ARITH == MUMPS_ARITH_z
498 if(im_rhs_sparse){
499 for(i=0;i<dmumps_par->nz_rhs;i++){
500 ((dmumps_par->rhs_sparse)[i]).r = rhs_sparse[i];
501 ((dmumps_par->rhs_sparse)[i]).i = im_rhs_sparse[i];
502 }
503 }else{
504 for(i=0;i<dmumps_par->nz_rhs;i++){
505 ((dmumps_par->rhs_sparse)[i]).r = rhs_sparse[i];
506 ((dmumps_par->rhs_sparse)[i]).i = 0.0;
507 }
508 }
509#else
510 for(i=0;i<dmumps_par->nz_rhs;i++){
511 (dmumps_par->rhs_sparse)[i] = rhs_sparse[i];
512 }
513#endif
514 }
515 }
516
517 if(dmumps_par->size_schur > 0){
518 if (dofactorize) {
519 MYMALLOC((dmumps_par->schur),((dmumps_par->size_schur)*(dmumps_par->size_schur)),double2);
520 }
521 dmumps_par->icntl[18] = 1;
522 }else{
523 dmumps_par->icntl[18] = 0;
524 }
525 /* Reduced RHS */
526 if ( dmumps_par->size_schur > 0 && dosolve ) {
527 if ( dmumps_par->icntl[26-1] == 2 ) {
528 /* REDRHS is on input */
529 tmp_m= mxGetM(REDRHS_IN);
530 tmp_n= mxGetN(REDRHS_IN);
531 if (tmp_m != dmumps_par->size_schur || tmp_n != dmumps_par->nrhs) {
532 mexErrMsgTxt ("bad dimensions for REDRHS in mumpsmex.c");
533 }
534 EXTRACT_CMPLX_FROM_MATLAB_TOPTR(REDRHS_IN,(dmumps_par->redrhs),double,((int)tmp_m*tmp_n));
535 dmumps_par->lredrhs=dmumps_par->size_schur;
536 }
537 if ( dmumps_par->icntl[26-1] == 1 ) {
538 /* REDRHS on output. Must be allocated before the call */
539 MYFREE(dmumps_par->redrhs);
540 if(!(dmumps_par->redrhs=(double2 *)malloc((dmumps_par->size_schur*dmumps_par->nrhs)*sizeof(double2)))){
541 mexErrMsgTxt("malloc redrhs failed in intmumpsc.c");
542 }
543 }
544 }
545 dmumps_c(dmumps_par);
546 }
547 }
548 if(nlhs > 0){
549 EXTRACT_FROM_C_TO_MATLAB( INFO_OUT ,(dmumps_par->infog),80);
550 EXTRACT_FROM_C_TO_MATLAB( RINFO_OUT ,(dmumps_par->rinfog),40);
551 /* A-1 on output */
552 if ( dmumps_par->icntl[30-1] != 0 && dosolve ) {
553 RHS_OUT = mxCreateSparse(dmumps_par->n, dmumps_par->n,dmumps_par->nz_rhs,mxREAL2);
554
555 irhs_ptr = mxGetJc(RHS_OUT);
556 irhs_sparse = mxGetIr(RHS_OUT);
557 for(j=0;j<dmumps_par->nrhs+1;j++){
558 irhs_ptr[j] = (mwIndex) ((dmumps_par->irhs_ptr)[j]-1);
559 }
560 ptr_matlab = mxGetPr(RHS_OUT);
561#if MUMPS_ARITH == MUMPS_ARITH_z
562 ptri_matlab = mxGetPi(RHS_OUT);
563#endif
564 for(i=0;i<dmumps_par->nz_rhs;i++){
565#if MUMPS_ARITH == MUMPS_ARITH_z
566 /* complex arithmetic */
567 ptr_matlab[i] = (dmumps_par->rhs_sparse)[i].r;
568 ptri_matlab[i] = (dmumps_par->rhs_sparse)[i].i;
569#else
570 /* real arithmetic */
571 ptr_matlab[i] = (dmumps_par->rhs_sparse)[i];
572#endif
573 irhs_sparse[i] = (mwIndex)((dmumps_par->irhs_sparse)[i]-1);
574 }
575
576 }
577 else if(dmumps_par->rhs && dosolve){
578 /* nbrhs may not have been set (case of null space) */
579 nbrhs=dmumps_par->nrhs;
580 RHS_OUT = mxCreateDoubleMatrix (dmumps_par->n,dmumps_par->nrhs,mxREAL2);
581 ptr_matlab = mxGetPr (RHS_OUT);
582#if MUMPS_ARITH == MUMPS_ARITH_z
583 ptri_matlab = mxGetPi (RHS_OUT);
584 for(j=0;j<dmumps_par->nrhs;j++){
585 posrhs = j*(int)n;
586 for(i=0;i<dmumps_par->n;i++){
587 ptr_matlab[posrhs+i]= (dmumps_par->rhs)[posrhs+i].r;
588 ptri_matlab[posrhs+i]= (dmumps_par->rhs)[posrhs+i].i;
589 }
590 }
591#else
592 for(j=0;j<dmumps_par->nrhs;j++){
593 posrhs = j*dmumps_par->n;
594 for(i=0;i<dmumps_par->n;i++){
595 ptr_matlab[posrhs+i]= (dmumps_par->rhs)[posrhs+i];
596 }
597 }
598#endif
599 }else{
600 EXTRACT_CMPLX_FROM_C_TO_MATLAB( RHS_OUT,(dmumps_par->rhs),1);
601 }
602
603 ptr_int = (int *)dmumps_par;
604 inst_address = (size_t) ptr_int;
605 EXTRACT_FROM_C_TO_MATLAB( INST_OUT ,&inst_address,1);
606 EXTRACT_FROM_C_TO_MATLAB( PIVNUL_LIST,dmumps_par->pivnul_list,dmumps_par->infog[27]);
607 EXTRACT_FROM_C_TO_MATLAB( PERM_OUT ,dmumps_par->sym_perm,dmumps_par->n);
608 EXTRACT_FROM_C_TO_MATLAB( UNS_PERM ,dmumps_par->uns_perm,dmumps_par->n);
609 EXTRACT_FROM_C_TO_MATLAB( ICNTL_OUT ,dmumps_par->icntl,60);
610 EXTRACT_FROM_C_TO_MATLAB( CNTL_OUT ,dmumps_par->cntl,15);
611 EXTRACT_FROM_C_TO_MATLAB( ROWSCA_OUT ,dmumps_par->rowsca,dmumps_par->n);
612 EXTRACT_FROM_C_TO_MATLAB( COLSCA_OUT ,dmumps_par->colsca,dmumps_par->n);
613 EXTRACT_FROM_C_TO_MATLAB( KEEP_OUT ,dmumps_par->keep,500);
614 EXTRACT_FROM_C_TO_MATLAB( DKEEP_OUT ,dmumps_par->dkeep,230);
615
616 if(dmumps_par->size_schur > 0 && dofactorize){
617 SCHUR_OUT = mxCreateDoubleMatrix(dmumps_par->size_schur,dmumps_par->size_schur,mxREAL2);
618 ptr_matlab = mxGetPr (SCHUR_OUT);
619#if MUMPS_ARITH == MUMPS_ARITH_z
620 ptri_matlab = mxGetPi (SCHUR_OUT);
621 for(i=0;i<dmumps_par->size_schur;i++){
622 pos = i*(dmumps_par->size_schur);
623 for(j=0;j<dmumps_par->size_schur;j++){
624 ptr_matlab[j+pos] = ((dmumps_par->schur)[j+pos]).r;
625 ptri_matlab[j+pos] = ((dmumps_par->schur)[j+pos]).i;
626 }
627 }
628#else
629 for(i=0;i<dmumps_par->size_schur;i++){
630 pos = i*(dmumps_par->size_schur);
631 for(j=0;j<dmumps_par->size_schur;j++){
632 ptr_matlab[j+pos] = (dmumps_par->schur)[j+pos];
633 }
634 }
635#endif
636 }else{
637 SCHUR_OUT = mxCreateDoubleMatrix(1,1,mxREAL2);
638 ptr_matlab = mxGetPr (SCHUR_OUT);
639 ptr_matlab[0] = -9999;
640#if MUMPS_ARITH == MUMPS_ARITH_z
641 ptr_matlab = mxGetPi (SCHUR_OUT);
642 ptr_matlab[0] = -9999;
643#endif
644 }
645 /* REDRHS on output */
646 if ( dmumps_par->icntl[26-1]==1 && dmumps_par->size_schur > 0 && dosolve ) {
647 REDRHS_OUT = mxCreateDoubleMatrix(dmumps_par->size_schur,dmumps_par->nrhs,mxREAL2);
648 ptr_matlab = mxGetPr(REDRHS_OUT);
649#if MUMPS_ARITH == MUMPS_ARITH_z
650 ptri_matlab = mxGetPi (REDRHS_OUT);
651#endif
652 for(i=0;i<dmumps_par->nrhs*dmumps_par->size_schur;i++){
653#if MUMPS_ARITH == MUMPS_ARITH_z
654 ptr_matlab[i] = ((dmumps_par->redrhs)[i]).r;
655 ptri_matlab[i] = ((dmumps_par->redrhs)[i]).i;
656#else
657 ptr_matlab[i] = ((dmumps_par->redrhs)[i]);
658#endif
659 }
660 }else{
661 REDRHS_OUT = mxCreateDoubleMatrix(1,1,mxREAL2);
662 ptr_matlab = mxGetPr (REDRHS_OUT);
663 ptr_matlab[0] = -9999;
664#if MUMPS_ARITH == MUMPS_ARITH_z
665 ptr_matlab = mxGetPi (REDRHS_OUT);
666 ptr_matlab[0] = -9999;
667#endif
668 }
669
670 MYFREE(dmumps_par->redrhs);
671 MYFREE(dmumps_par->schur);
672 MYFREE(dmumps_par->irhs_ptr);
673 MYFREE(dmumps_par->irhs_sparse);
674 MYFREE(dmumps_par->rhs_sparse);
675 MYFREE(dmumps_par->rhs);
676 }
677}
dmumps_c
void MUMPS_CALL dmumps_c(DMUMPS_STRUC_C *dmumps_par)
ROWSCA_IN
#define ROWSCA_IN
Definition mumpsmex.c:38
EXTRACT_SCALING_FROM_MATLAB_TOPTR
#define EXTRACT_SCALING_FROM_MATLAB_TOPTR(mxcomponent, mumpspointer, is_a_pointer_from_mumps, length)
Definition mumpsmex.c:91
RINFO_OUT
#define RINFO_OUT
Definition mumpsmex.c:48
RHS_OUT
#define RHS_OUT
Definition mumpsmex.c:49
REDRHS_IN
#define REDRHS_IN
Definition mumpsmex.c:42
CNTL_IN
#define CNTL_IN
Definition mumpsmex.c:35
INST_OUT
#define INST_OUT
Definition mumpsmex.c:50
INFO_OUT
#define INFO_OUT
Definition mumpsmex.c:47
ICNTL_IN
#define ICNTL_IN
Definition mumpsmex.c:34
A_IN
#define A_IN
Definition mumpsmex.c:45
PERM_OUT
#define PERM_OUT
Definition mumpsmex.c:54
DKEEP_IN
#define DKEEP_IN
Definition mumpsmex.c:44
EXTRACT_FROM_MATLAB_TOVAL
#define EXTRACT_FROM_MATLAB_TOVAL(mxcomponent, mumpsvalue)
Definition mumpsmex.c:111
PERM_IN
#define PERM_IN
Definition mumpsmex.c:36
ROWSCA_OUT
#define ROWSCA_OUT
Definition mumpsmex.c:59
DKEEP_OUT
#define DKEEP_OUT
Definition mumpsmex.c:61
RHS_IN
#define RHS_IN
Definition mumpsmex.c:39
KEEP_IN
#define KEEP_IN
Definition mumpsmex.c:43
DMUMPS_free
void DMUMPS_free(DMUMPS_STRUC_C **dmumps_par)
Definition mumpsmex.c:174
VAR_SCHUR
#define VAR_SCHUR
Definition mumpsmex.c:40
JOB
#define JOB
Definition mumpsmex.c:33
COLSCA_IN
#define COLSCA_IN
Definition mumpsmex.c:37
PIVNUL_LIST
#define PIVNUL_LIST
Definition mumpsmex.c:53
CNTL_OUT
#define CNTL_OUT
Definition mumpsmex.c:57
UNS_PERM
#define UNS_PERM
Definition mumpsmex.c:55
REDRHS_OUT
#define REDRHS_OUT
Definition mumpsmex.c:52
INST
#define INST
Definition mumpsmex.c:41
SYM
#define SYM
Definition mumpsmex.c:32
ICNTL_OUT
#define ICNTL_OUT
Definition mumpsmex.c:56
COLSCA_OUT
#define COLSCA_OUT
Definition mumpsmex.c:58
SCHUR_OUT
#define SCHUR_OUT
Definition mumpsmex.c:51
EXTRACT_FROM_C_TO_MATLAB
#define EXTRACT_FROM_C_TO_MATLAB(mxcomponent, mumpspointer, length)
Definition mumpsmex.c:117
DMUMPS_alloc
void DMUMPS_alloc(DMUMPS_STRUC_C **dmumps_par)
Definition mumpsmex.c:204
KEEP_OUT
#define KEEP_OUT
Definition mumpsmex.c:60
EXTRACT_FROM_MATLAB_TOARR
#define EXTRACT_FROM_MATLAB_TOARR(mxcomponent, mumpsarray, type, length)
Definition mumpsmex.c:101
EXTRACT_FROM_MATLAB_TOPTR
#define EXTRACT_FROM_MATLAB_TOPTR(mxcomponent, mumpspointer, type, length)
Definition mumpsmex.c:74
r
r
Definition schur_example.m:55
n
n
Definition schur_example.m:9
DMUMPS_STRUC_C::icntl
MUMPS_INT icntl[60]
Definition dmumps_c.h:46
DMUMPS_STRUC_C::rhs
DMUMPS_COMPLEX * rhs
Definition dmumps_c.h:95
DMUMPS_STRUC_C::colsca
DMUMPS_REAL * colsca
Definition dmumps_c.h:89
DMUMPS_STRUC_C::dkeep
DMUMPS_REAL dkeep[230]
Definition dmumps_c.h:49
DMUMPS_STRUC_C::infog
MUMPS_INT infog[80]
Definition dmumps_c.h:100
DMUMPS_STRUC_C::redrhs
DMUMPS_COMPLEX * redrhs
Definition dmumps_c.h:95
DMUMPS_STRUC_C::nz
MUMPS_INT nz
Definition dmumps_c.h:58
DMUMPS_STRUC_C::job
MUMPS_INT job
Definition dmumps_c.h:44
DMUMPS_STRUC_C::irn
MUMPS_INT * irn
Definition dmumps_c.h:60
DMUMPS_STRUC_C::n
MUMPS_INT n
Definition dmumps_c.h:51
DMUMPS_STRUC_C::a
DMUMPS_COMPLEX * a
Definition dmumps_c.h:62
DMUMPS_STRUC_C::sym
MUMPS_INT sym
Definition dmumps_c.h:44
DMUMPS_STRUC_C::schur
DMUMPS_COMPLEX * schur
Definition dmumps_c.h:111
DMUMPS_STRUC_C::rinfog
DMUMPS_REAL rinfog[40]
Definition dmumps_c.h:101
DMUMPS_STRUC_C::rhs_sparse
DMUMPS_COMPLEX * rhs_sparse
Definition dmumps_c.h:95
DMUMPS_STRUC_C::colsca_from_mumps
MUMPS_INT colsca_from_mumps
Definition dmumps_c.h:91
DMUMPS_STRUC_C::rowsca_from_mumps
MUMPS_INT rowsca_from_mumps
Definition dmumps_c.h:92
DMUMPS_STRUC_C::perm_in
MUMPS_INT * perm_in
Definition dmumps_c.h:82
DMUMPS_STRUC_C::cntl
DMUMPS_REAL cntl[15]
Definition dmumps_c.h:48
DMUMPS_STRUC_C::size_schur
MUMPS_INT size_schur
Definition dmumps_c.h:109
DMUMPS_STRUC_C::uns_perm
MUMPS_INT * uns_perm
Definition dmumps_c.h:86
DMUMPS_STRUC_C::irhs_sparse
MUMPS_INT * irhs_sparse
Definition dmumps_c.h:96
DMUMPS_STRUC_C::nz_alloc
MUMPS_INT nz_alloc
Definition dmumps_c.h:54
DMUMPS_STRUC_C::rowsca
DMUMPS_REAL * rowsca
Definition dmumps_c.h:90
DMUMPS_STRUC_C::keep
MUMPS_INT keep[500]
Definition dmumps_c.h:47
DMUMPS_STRUC_C::listvar_schur
MUMPS_INT * listvar_schur
Definition dmumps_c.h:110
DMUMPS_STRUC_C::nz_rhs
MUMPS_INT nz_rhs
Definition dmumps_c.h:97
DMUMPS_STRUC_C::irhs_ptr
MUMPS_INT * irhs_ptr
Definition dmumps_c.h:96
DMUMPS_STRUC_C::lrhs
MUMPS_INT lrhs
Definition dmumps_c.h:97
DMUMPS_STRUC_C::par
MUMPS_INT par
Definition dmumps_c.h:44
DMUMPS_STRUC_C::pivnul_list
MUMPS_INT * pivnul_list
Definition dmumps_c.h:105
DMUMPS_STRUC_C::jcn
MUMPS_INT * jcn
Definition dmumps_c.h:61
DMUMPS_STRUC_C::sym_perm
MUMPS_INT * sym_perm
Definition dmumps_c.h:85
DMUMPS_STRUC_C::nrhs
MUMPS_INT nrhs
Definition dmumps_c.h:97
DMUMPS_STRUC_C::lredrhs
MUMPS_INT lredrhs
Definition dmumps_c.h:97