SPARSESUB
ModelingSolves the linear system of equations of each Newton iteration. The linear system is defined as $\mathit{J}\u2206q=F(q,\dot{q},\ddot{q},t)$ , where $\mathit{J}$ is a sparse, unsymmetrical Jacobian matrix, $\u2206q$ is the difference in states, and $F(q,\dot{q},\ddot{q},t)$ is the righthand side.
Use
An example of the corresponding Param_Simulation
statement:
<Param_Simulation
Linsolver = "USER"
/>
<Param_Sparse
usrsub_param_string = "USER([par1,par2,…,parN])"
usrsub_dll_name = "mySub"
usrsub_fnc_name = "SPARSESUB"
/>
Format
 Fortran Calling Syntax

SUBROUTINE SPARSESUB (PAR, NPAR, N, NNZ, AX, AI, AP, JOBFLG, IFLAG, RHS, ERRFLAG)
 C/C++ Calling Syntax

void STDCALL SPARSESUB (double *PAR, int *NPAR, int *N, int *NNZ, double *AX, int *AI, int *AP, int *JOBFLG, int *IFLAG, double *RHS, int *ERRFLAG)
Attributes
 PAR
 [double precision]
 NPAR
 [integer]
 N
 [integer]
 NNZ
 [integer]
 AX
 [double]
 AI
 [integer]
 AP
 [integer]
 JOBFLG
 [integer]
 IFLAG
 [integer]
 RHS
 [double]
Output
 ERRFLG
 [integer]
Example
The following C++ user subroutine integrates the PARDISO
solver.
void STDCALL SPARSESUB(double *PAR, int *NPAR, int *N, int *NNZ, double *AX, int *AI, int *AP, int *JOBFLG, int *IFLAG, double *RHS, int *ERRFLAG)
{
if (*IFLAG == 1)
{
_iwork = new int[64]; // Use as IPARM
memset(_iwork, 0, 64 * sizeof(int));
_x = new double[*N]; // Use as X
memset(_x, 0, *N * sizeof(double));
for (int i = 0; i < 64; i++)
{
_jac_pt[i] = 0;
}
_iwork[0] = 1; /* No solver default */
_iwork[1] = 2; /* Fillin reordering from METIS */
_iwork[2] = 1; /* Numbers of processors, value of OMP_NUM_THREADS */
_iwork[3] = 0; /* No iterativedirect algorithm */
_iwork[4] = 0; /* No user fillin reducing permutation */
_iwork[5] = 0; /* Write solution into x */
_iwork[6] = 0; /* Not in use */
_iwork[7] = 2; /* Max numbers of iterative refinement steps */
_iwork[8] = 0; /* Not in use */
_iwork[9] = 13; /* Perturb the pivot elements with 1E13 */
_iwork[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
_iwork[11] = 0; /* Conjugate transposed/transpose solve */
_iwork[12] = 1; /* Maximum weighted matching algorithm is switchedon (default for nonsymmetric) */
_iwork[13] = 0; /* Output: Number of perturbed pivots */
_iwork[14] = 0; /* Not in use */
_iwork[15] = 0; /* Not in use */
_iwork[16] = 0; /* Not in use */
_iwork[17] = 1; /* Output: Number of nonzeros in the factor LU */
_iwork[18] = 1; /* Output: Mflops for LU factorization */
_iwork[19] = 0; /* Output: Numbers of CG Iterations */
// _iwork[20] = 3; /* CG/CGS diagnostics */
_iwork[34] = 1; /* Use zerobased indexing: columns and rows indexing in arrays ia , ja , and perm starts from 0 (Cstyle indexing)*/
}
else
{
int neqs = *N;
MKL_INT mtype = 11; /* Real unsymmetric matrix */
MKL_INT nrhs = 1; /* Number of right hand sides. */
MKL_INT maxfct = 1; /* Maximum number of numerical factorizations. */
MKL_INT mnum = 1; /* Which factorization to use. */
MKL_INT phase = 22;
MKL_INT error = 0; /* Initialize error flag */
MKL_INT msglvl = 0; /* Print statistical information in file */
MKL_INT idum; /* Integer dummy. */
if (*JOBFLG == 1)
phase = 11; // Symbolic Factorization
else if (*JOBFLG == 2  *JOBFLG == 3)
phase = 22; // Numerical Factorization/ReFactorization
else if (*JOBFLG == 4)
phase = 33; // Solve RHS
PARDISO(_jac_pt, &maxfct, &mnum, &mtype, &phase,
&neqs, AX, AI, AP, &idum, &nrhs, _iwork, &msglvl, RHS, _x, &error);
if (*JOBFLG == 4)
memcpy(RHS, _x, neqs * sizeof(double));
if (error != 0)
*ERRFLAG = 1;
else
*ERRFLAG = 0;
}
}