SPARSESUB

ModelingSolves the linear system of equations of each Newton iteration. The linear system is defined as J ∆ q = - F ( q , q ˙ , q ¨ , t ) , where J is a sparse, unsymmetrical Jacobian matrix, ∆ q is the difference in states, and - F ( q , q ˙ , q ¨ , t ) is the right-hand 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]
An array that contains the constant arguments from the list provided in the user-defined statement.
NPAR
[integer]
The number of entries in the PAR array.
N
[integer]
Number of equations in the sparse linear system of equations J ∆ q = - F ( q , q ˙ , q ¨ , t ) . The matrix J is of size N x N.
NNZ
[integer]
Number of nonzero entries in the sparse matrix.
AX
[double]
Nonzero values of the Jacobian matrix J . The array AX contains the values corresponding to the indices in AP. The size and order of AX is the same as that of AP. The coefficients can be either real or complex. The matrix must be stored in a compressed sparse row format with increasing values of AP for each row.
AI
[integer]
Integer array of size N+1. It points to the first column index of row i in the array AP in a compressed sparse row format.
Note: Row and column numbers start from 1.
AP
[integer]
Contains the column indices of the sparse matrix J stored in a compressed sparse row format. The indices in each row must be sorted in increasing order.
JOBFLG
[integer]
Defines the job to execute:
  • 1 = Symbolic factorization
  • 2 = Numerical factorization
  • 3 = Numerical refactorization
  • 4 = Solve RHS
IFLAG
[integer]
An integer variable that MotionSolve sets to 1 when the user subroutine is called during MotionSolve initialization. The IFLAG is set to zero at all other times.
RHS
[double]
Defines the right-hand side − F q , q ˙ , q ¨ , t MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqGHsislcaWGgbWaaeWaa8aabaWdbiaadghacaGGSaGabmyCa8aa gaGaa8qacaGGSaWdamaaxacabaWdbiaadghaaSWdaeqabaWdbiaacI kaaaGccaGGSaGaamiDaaGaayjkaiaawMcaaaaa@4144@ of the Newton iteration. It is a vector of size N.

Output

ERRFLG
[integer]
Information about the call status. A nonzero value indicates failure.

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;       /* Fill-in reordering from METIS */
        _iwork[2] = 1;       /* Numbers of processors, value of OMP_NUM_THREADS */
        _iwork[3] = 0;       /* No iterative-direct algorithm */
        _iwork[4] = 0;       /* No user fill-in 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 1E-13 */
        _iwork[10] = 1;      /* Use nonsymmetric permutation and scaling MPS */
        _iwork[11] = 0;      /* Conjugate transposed/transpose solve */
        _iwork[12] = 1;      /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
        _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 zero-based indexing: columns and rows indexing in arrays ia , ja , and perm starts from 0 (C-style 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;
    }
}