3. System Vectors and Matrices - SysVector, SysMatrix

The two modules SysVector and SysMatrix are designed to manage and operate upon global degree of freedom (dof) system vectors and matrices. System vectors result, for example, from the assembly of individual finite element vectors, such as consistent applied load vectors. The SysVector module provides functions for assembly, query and a variety of arithmetic operations. System matrices result from the assembly of finite element matrices, such as stiffness matrices. System matrices may be symmetric, non-symmetric or diagonal and held in sparse, full, or iterative data structures. Functions exist for assembly, direct and iterative solution, eigenvalue and eigenvector extraction, etc.

3.1. System Vectors - SysVector

The SysVector module manages and operates upon global dof system vectors. Examples of system vectors are assembled load vectors, displacement solution vectors, eigenvectors, residual force vectors, etc. A SysVector object may require a significant amount of memory to hold values at each global dof. The vector values are represented internally as a contiguous vector of double precision numbers. The length of the vector is equal to the number of global dof. SysVector objects are used as arguments to SysMatrix functions for performing system matrix-vector products, etc. The methods associated with a SysVector object are the following.

Instance a SysVector object using vfs_SysVectorBegin(). Once a SysVector object is instanced, specify whether the system refers to a complex matrix, if needed, using vfs_SysVectorPre(). Set the number of global dof to be managed using vfs_SysVectorDef(). Initially all global dof values are set to 0.

All values in a SysVector object may be initialized to some user defined real value using vfs_SysVectorInit(). Values may be set into a SysVector object at selected degrees of freedom using vfs_SysVectorScatter(). Values may be added into a SysVector object at selected degrees of freedom using vfs_SysVectorAssem(). Note that the values should be initialized to zero before using vfs_SysVectorAssem().

Random values may be generated using vfs_SysVectorRandom().

Various operations can be performed upon system vectors such as dot product addition, scaling, etc using vfs_SysVectorDot(), vfs_SysVectorAdd(), vfs_SysVectorScale(), etc.

Values may be selectively retrieved from a system vector at selected global dof numbers using vfs_SysVectorGather().

SysVector objects are entered as arguments to a SysMatrix object for performing matrix-vector products, back substitution on factored system matrices, etc. Destroy an instance of a SysVector object using vfs_SysVectorEnd().

If Gather/Scatter/Assem operations are performed in parallel over several threads then the complex mode can be specified on each thread using vfs_SysVectorSetThreadComplexMode().

3.2. Function Descriptions

The currently available SysVector functions are described in detail in this section.

vfs_SysVector *vfs_SysVectorBegin(void)

create an instance of a SysVector object

Create an instance of a SysVector object. Memory is allocated for the object private data and the pointer to the object is returned.

Destroy an instance of a SysVector object using 

void vfs_SysVectorEnd (vfs_SysVector *sysvector)

Return the current value of a SysVector object error flag using 

Vint vfs_SysVectorError (vfs_SysVector *sysvector)

Make a copy of a SysVector object. The private data from the fromsysvector object is copied to the sysvector object. Any previous private data in sysvector is lost. 

void vfs_SysVectorCopy (vfs_SysVector *sysvector,
                        vfs_SysVector *fromsysvector)

Returns:

The function returns a pointer to the newly created SysVector object. If the object creation fails, NULL is returned.

void vfs_SysVectorEnd(vfs_SysVector *p)

destroy an instance of a SysVector object

See vfs_SysVectorBegin()

Vint vfs_SysVectorError(vfs_SysVector *p)

return the current value of a SysVector object error flag

See vfs_SysVectorBegin()

void vfs_SysVectorDef(vfs_SysVector *p, Vint neq)

define number of global dof

Specify the number of global dof. The vector values are undefined.

Inquire of a defined neq as an output argument using 

void vfe_SysVectorInq (vfe_SysVector *sysvector,
                       Vint *neq)

Errors

  • SYS_ERROR_VALUE is generated if neq is less than or equal to zero.

  • SYS_ERROR_MEMORY is generated if a memory allocation failure occurs.

Parameters:

  • p – Pointer to SysVector object.

  • neq – Number of global dof

void vfs_SysVectorInq(vfs_SysVector *p, Vint *neq)

inquire of a defined neq as an output argument

See vfs_SysVectorDef()

void vfs_SysVectorInit(vfs_SysVector *p, Vdouble v)

initialize vector to a given value

Initialize all entries in the system vector to v.

Parameters:

  • p – Pointer to SysVector object.

  • v – Initialization value

void vfs_SysVectorRandom(vfs_SysVector *p)

generate a random vector

Initialize a vector to random values in the interval [-1.,1.]. If several random vectors are to be generated they should be generated in a single SysVector object and moved to another SysVector object. If each random vector were to be generated in separate SysVector objects they would all be random but identical.

Parameters:

p – Pointer to SysVector object.

void vfs_SysVectorAssem(vfs_SysVector *p, Vint nedofs, Vint lm[], Vdouble v[])

assemble an element vector

Assemble the element vector v. The values in v are added to the existing global dof values. The number of element degrees of freedom for the element is nedofs. The global dof numbers are input in the lm vector.

Errors

SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

Parameters:

  • p – Pointer to SysVector object.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • v – Element vector

void vfs_SysVectorThreadAssem(vfs_SysVector *p, Vint ithread, Vint nedofs, Vint lm[], Vdouble v[])

assemble an element vector in a thread

Assemble the element vector v. The values in v are added to the existing global dof values. The number of element degrees of freedom for the element is nedofs. The global dof numbers are input in the lm vector.

If a complex mode is specified with vfs_SysVectorSetThreadComplexMode() then the array v refers to the mode specified.

Errors

SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ithread – Thread number starting from 1.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • v – Element vector

void vfs_SysVectorGather(vfs_SysVector *p, Vint ndofs, Vint lm[], Vdouble v[])

gather values from vector

Gather ndofs values at global dof numbers lm and return in vector v.

Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers

  • v – [out] Gathered values

void vfs_SysVectorThreadGather(vfs_SysVector *p, Vint ithread, Vint ndofs, Vint lm[], Vdouble v[])

gather values from vector in thread

Gather ndofs values at global dof numbers lm and return in vector v.

If a complex mode is specified with vfs_SysVectorSetThreadComplexMode() then the array v refers to the mode specified.

Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ithread – Thread number starting from 1.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers

  • v – [out] Gathered values

void vfs_SysVectorScatter(vfs_SysVector *p, Vint ndofs, Vint lm[], Vdouble v[])

scatter values to vector

Scatter ndofs values in vector v at global dof numbers lm to the system vector sysvector. In matrix form, given 

sysvector = [ a_1 a_2 ... a_N ]
then 
a_lm[0] = v[0]
a_lm[1] = v[1]
     ...
a_lm[ndofs-1] = v[ndofs-1]
Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers

  • v – Scattered values

void vfs_SysVectorThreadScatter(vfs_SysVector *p, Vint ithread, Vint ndofs, Vint lm[], Vdouble v[])

scatter values to vector in thread

Scatter ndofs values in vector v at global dof numbers lm to the system vector sysvector.

If a complex mode is specified with vfs_SysVectorSetThreadComplexMode() then the array v refers to the mode specified.

Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ithread – Thread number starting from 1.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers

  • v – Scattered values

void vfs_SysVectorScatterVal(vfs_SysVector *p, Vint ndofs, Vint lm[], Vdouble v)

scatter a constant value to vector

Scatter the value v to the ndofs values into the system vector sysvector at global dof numbers lm.

Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysVector object.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers

  • v – Constant scattered value

void vfs_SysVectorDot(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vdouble *s)

compute dot product of two vectors

Compute the dot product of sysvector and sysvectorx. A single SysVector object can be safely used as more than one argument.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system vector is not equal.

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object.

  • s – [out] Dot product

void vfs_SysVectorScale(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vdouble s)

scale a vector by a value

Scale the values of vector sysvectorx and set in sysvector. The operation replaces any previous values in sysvector with the scaled values of sysvectorx. A single SysVector object can be safely used as more than one argument.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system vector is not equal.

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object.

  • s – Scale factor

void vfs_SysVectorNorm(vfs_SysVector *p, Vint type, Vdouble *value)

compute the norm of a vector

Compute and return in value the vector norm specified by type. Use vfs_SysVectorLInfNorm() to compute the L-Infinity norm if the global degree of freedom number at which the maximum component absolute value occurs is required.

Parameters:

  • p – Pointer to SysVector object.

  • type – Type of vector norm to compute 

    x=SYSVECTOR_NORM_L1         Sum of absolute values
 =SYSVECTOR_NORM_L2         Square root of sum of squares
 =SYSVECTOR_NORM_LINF       Maximum component absolute value

  • value – [out] Norm value

void vfs_SysVectorLInfNorm(vfs_SysVector *p, Vdouble *value, Vint *dof)

compute the L-Infinity norm of a vector

Compute and return in value the L-infinity norm. The global degree of freedom number at which the maximum component absolute value occurs is returned in dof.

Parameters:

  • p – Pointer to SysVector object.

  • value – [out] L-Infinity norm value

  • dof – [out] Global dof number

void vfs_SysVectorAdd(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vdouble s, vfs_SysVector *sysvectory)

add two vectors

Add two vectors, sysvectorx and sysvectory to sysvector . The operation replaces any previous values in SysVector with the sum of sysvectorx and s times sysvectory. A single SysVector object can be safely used as more than one argument.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system vector is not equal.

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object.

  • s – Scale factor

  • sysvectory – Pointer to SysVector object.

void vfs_SysVectorReduce(vfs_SysVector *p, Vint nvec, vfs_SysVector *sysvectort[], vfs_SysVector *sysvectorq)

reduce a system vector

Reduce the system vector, sysvector, by nvec vectors, sysvectort[], to yield a reduced system vector, sysvectorq. The order of the reduced vector, sysvectorq, must be nvec. In matrix form, given 

sysvectorq = [ a_1 a_2 ... a_N ]
then 
a_1 = sysvector * sysvectort[0]
a_2 = sysvector * sysvectort[1]
        ...
a_N = sysvector * sysvectort[N-1]
Errors

  • SYS_ERROR_VALUE is generated if an improper nvec is specified.

  • SYS_ERROR_OPERATION is generated if the number of global dof in sysvectort and sysvector are not equal

Parameters:

  • p – Pointer to SysVector object.

  • nvec – Number of system vectors 0 < nvec

  • sysvectort – Poinsysvectorter to array of nvec SysVector objects

  • sysvectorq – Pointer to reduced SysVector object.

void vfs_SysVectorExpand(vfs_SysVector *p, Vint nvec, vfs_SysVector *sysvectort[], vfs_SysVector *sysvectorq)

expand a system vector

Expand the system vector, sysvector, by nvec vectors, sysvectort[], to yield an expanded system vector, sysvectorq. The order of the vector to be expanded, sysvector, must be nvec. In matrix form, given 

sysvector = [ a_1 a_2 ... a_N ]
then 
sysvectorq = [ a_1 * sysvectort[0] + a_2 * sysvectort[1] + ... + a_N * sysvectort[N-1]
Errors

  • SYS_ERROR_VALUE is generated if an improper nvec is specified.

  • SYS_ERROR_OPERATION is generated if the number of global dof in sysvectort and sysvectorq are not equal

Parameters:

  • p – Pointer to SysVector object.

  • nvec – Number of system vectors 0 < nvec

  • sysvectort – Pointer to array of nvec SysVector objects

  • sysvectorq – Pointer to expanded SysVector object.

void vfs_SysVectorCopy(vfs_SysVector *p, vfs_SysVector *q)

make a copy of a SysVector object

See vfs_SysVectorBegin()

void vfs_SysVectorMove(vfs_SysVector *p, vfs_SysVector *sysvectorx)

move vector values to another vector

Move the values in sysvectorx to sysvector. This operation is different from the deep copy operation vfs_SysVectorCopy() which copies the vector values as well as the number of global dof.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system vector is not equal.

Parameters:
void vfs_SysVectorProd(vfs_SysVector *p, vfs_SysVector *sysvectorx, vfs_SysVector *sysvectory)

compute the product of two vectors

Compute the product of two vectors, sysvectorx and sysvectory, component-wise and place the result in sysvector. Use vfs_SysVectorAddProd() to add the product to sysvector.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in sysvectorx, sysvectory and sysvector are not equal

Parameters:
void vfs_SysVectorAddProd(vfs_SysVector *p, Vdouble s, vfs_SysVector *sysvectorx, vfs_SysVector *sysvectory)

compute and add the product of two vectors

Compute the product of two vectors, sysvectorx and sysvectory, component-wise scaled by s and add the result to sysvector.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in sysvectorx, sysvectory and sysvector are not equal

Parameters:

  • p – Pointer to SysVector object.

  • s – Scale factor

  • sysvectorx – Pointer to SysVector object.

  • sysvectory – Pointer to SysVector object.

void vfs_SysVectorInverse(vfs_SysVector *p, vfs_SysVector *sysvectorx)

compute the inverse of a vector

Compute the inverse of vector sysvectorx component-wise and place in sysvector.

Errors

  • SYS_ERROR_OPERATION is generated if the number of global dof in each system vector is not equal.

  • SYS_ERROR_COMPUTE is generated if sysvectorx has a zero component.

Parameters:
void vfs_SysVectorSquareRoot(vfs_SysVector *p, vfs_SysVector *sysvectorx)

computes square root of vector

Use vfs_SysVectorSquareRoot() to take the square root of every component in sysvectorx, storing the result in sysvector.

Errors

SYS_ERROR_COMPUTE is generated if any component in sysvectorx is negative.

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object to take square root of.

void vfs_SysVectorMaskGather(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vint n, Vint mask[])

gather values from vector with mask

Gather into the first n entries of system vector sysvector the components of sysvectorx referenced in array mask. In matrix form, given 

sysvector  = [ a_1 a_2 ... a_N ]
sysvectorx = [ x_1 x_2 ... x_N ]
then 
a_1 = x_mask[0]
a_2 = x_mask[1]
      ...
a_[n-1] = x_mask[n-1]

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object where data will be copied from.

  • n – Number of components to copy

  • mask – Vector of 0-based components to use when referencing system vector sysvectorx.

void vfs_SysVectorMaskScatter(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vint n, Vint mask[])

scatter values to vector with mask

Gather the first n entries of system vector sysvectorx into the components of sysvector referenced in array mask. In matrix form, given 

sysvector  = [ a_1 a_2 ... a_N ]
sysvectorx = [ x_1 x_2 ... x_N ]
then 
a_mask[0] = x_1
a_mask[1] = x_2
      ...
a_mask[n-1] = x_[n-1]

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object where data will be copied from.

  • n – Number of components to copy

  • mask – Vector of 0-based components to use when referencing system vector sysvectorx.

void vfs_SysVectorMaskAdd(vfs_SysVector *p, vfs_SysVector *sysvectorx, Vdouble s, vfs_SysVector *sysvectory, Vint n, Vint mask[])

add two vectors with mask

Adds the first n entries of system vector sysvectorx to the scaled components of sysvectory referenced in array mask. In matrix form, given 

sysvector  = [ a_1 a_2 ... a_N ]
sysvectorx = [ x_1 x_2 ... x_N ]
sysvectory = [ y_1 y_2 ... y_n ]
then 
a_1 = x_1 + s * y_mask[0]
a_2 = x_2 + s * y_mask[1]
      ...
a_[n-1] = x_[n-1] + s * y_mask[n-1]

Parameters:

  • p – Pointer to SysVector object.

  • sysvectorx – Pointer to SysVector object containing first vector in sum.

  • s – Scaling parameter.

  • sysvectory – Pointer to SysVector object containing second vector in sum, whose components given by the 0-based mask array will be scaled by s during the sum.

  • n – Number of components to copy

  • mask – Vector of 0-based components to use when referencing system vector sysvectory.

void vfs_SysVectorPre(vfs_SysVector *p, Vint complextype)

specify real or complex vector

Use vfs_SysVectorPre() to specify whether the matrix is real or complex. For real vectors complextype takes the value SYS_DOUBLE. For complex vectors, complextype takes the value SYS_DOUBLECOMPLEX.

This call must be made before vfs_SysVectorDef() is called. Defaults to real.

Errors

SYS_ERROR_ENUM is issued if complextype is invalid.

Parameters:

  • p – Pointer to SysVector object.

  • complextype – Enumerated type for real or complex matrix 

    x=SYS_DOUBLE        Real matrix
 =SYS_DOUBLECOMPLEX Complex matrix

void vfs_SysVectorSetComplexMode(vfs_SysVector *p, Vint complexmode)

set complex mode

Specify floating point data mode.

The methods vfs_SysVectorAssem(), vfs_SysVectorGather(), and vfs_SysVectorScatter() all take arrays of floating point numbers. The values entered or returned can be declared as SYS_COMPLEX_REAL, SYS_COMPLEX_IMAGINARY, or SYS_COMPLEX_REALIMAGINARY. The data in the argument then takes the specified mode.

This is useful if only real or only imaginary data is being processed. The mode SYS_COMPLEX_REAL indicates that only one value, the real part, is being entered for each a_i in the array. In this case, the imaginary entries in the matrix remain unchanged. Likewise, the mode SYS_COMPLEX_IMAGINARY only affects the imaginary part. If mode is set to SYS_COMPLEX_REALIMAGINARY then two consecutive numbers in the array stand for the real and imaginary parts of the data. Defaults to SYS_COMPLEX_REAL.

Errors

  • SYS_ERROR_ENUM is generated if an improper complexmode is specified.

  • SYS_ERROR_OPERATION is generated if the matrix is not defined as complex.

Parameters:

  • p – Pointer to SysVector object.

  • complexmode – The type of data being passed through the Vdouble[] argument list. 

    x=SYS_COMPLEX_REAL          Data is real
 =SYS_COMPLEX_IMAGINARY     Data is imaginary
 =SYS_COMPLEX_REALIMAGINARY Data is complex

void vfs_SysVectorSetThreadComplexMode(vfs_SysVector *p, Vint ithread, Vint complexmode)

set complex mode in thread

Specify floating point data mode.

The methods vfs_SysVectorAssem(), vfs_SysVectorGather(), and vfs_SysVectorScatter() all take arrays of floating point numbers. The values entered or returned can be declared as SYS_COMPLEX_REAL, SYS_COMPLEX_IMAGINARY, or SYS_COMPLEX_REALIMAGINARY. The data in the argument then takes the specified mode. Each thread in parallel operations can take different values for the complexmode value.

This is useful if only real or only imaginary data is being processed. The mode SYS_COMPLEX_REAL indicates that only one value, the real part, is being entered for each a_i in the array. In this case, the imaginary entries in the matrix remain unchanged. Likewise, the mode SYS_COMPLEX_IMAGINARY only affects the imaginary part. If mode is set to SYS_COMPLEX_REALIMAGINARY then two consecutive numbers in the array stand for the real and imaginary parts of the data. Defaults to SYS_COMPLEX_REAL.

Errors

  • SYS_ERROR_ENUM is generated if an improper complexmode is specified.

  • SYS_ERROR_OPERATION is generated if the matrix is not defined as complex.

Parameters:

  • p – Pointer to SysVector object.

  • ithread – Thread number starting from 1.

  • complexmode – The type of data being passed through the Vdouble[] argument list. 

    x=SYS_COMPLEX_REAL          Data is real
 =SYS_COMPLEX_IMAGINARY     Data is imaginary
 =SYS_COMPLEX_REALIMAGINARY Data is complex

void vfs_SysVectorSetNumThreads(vfs_SysVector *p, Vint num)

set number of threads

Set the number of threads used for parallel execution of SysVector functions. If num is 1, then SysVector functions are run serially. If num is greater than or equal to 2 then num threads are used for parallel execution.

Errors

SYS_ERROR_VALUE is generated if num is less than 1.

Parameters:

  • p – Pointer to SysVector object.

  • num – Number of threads 

    x=1                   Default, serial execution.
 >1                   Number of threads for parallel execution

3.3. System Matrices - SysMatrix

The SysMatrix module manages and operates upon global dof system matrices. System matrices result, for example, from the assembly of finite element stiffness, mass or damping matrices. The SysMatrix module provides functions for assembly, direct and iterative solution, eigenvalue and eigenvector extraction, etc. System matrices may be symmetric, non-symmetric or diagonal and held in sparse, diagonal or full data structures. The methods associated with a SysMatrix object are the following.

Instance a SysMatrix object using vfs_SysMatrixBegin(). Once a SysMatrix object is instanced, specify whether the system refers to a complex matrix, if needed, using vfs_SysMatrixPre(). Set the number of global dof to be managed and the matrix type using vfs_SysMatrixDef(). The matrix type may be either sparse symmetric or non-symmetric, full symmetric or non-symmetric, diagonal, or structured for symmetric or non-symmetric iterative solution. Use full symmetric for small full matrices which result from some type of substructuring or projection process.

The SysMatrix module utilizes the BLAS routines for performance reasons. Use VKI_LIBAPI_BLASMKL_SEQUENTIAL or VKI_LIBAPI_BLASMKL_THREAD when compiling with Intel’s MKL BLAS routines for optimal performance. Use VKI_LIBAPI_BLASMKL_ILP64 if your application uses the 64-bit integer API for the MKL BLAS library. Use VKI_LIBAPI_OPENBLAS if the OpenBLAS are to be used.

See section Compiling and Linking a VfsTools Application. The complete, external BLAS functions are required for the Block Lanczos and AMLS eigensolvers.

The direct sparse solver used in SysMatrix may be changed from the default native solver. Two external solver packages are supported, Intel Pardiso and MUMPS. See section Compiling and Linking a VfsTools Application.

The assembled matrix is established in a SysMatrix object by assembly of individual finite element matrices and optional multipoint constraints. Multipoint constraints in general may be enforced using a variety of methods such as penalty functions, Lagrange multipliers or direct matrix transformation. If multipoint constraints are enforced using penalty functions or Lagrange multipliers then they are assembled as symmetric “element” matrices in the conventional way. Note that the direct sparse solvers currently do not support full pivoting which is necessary for robust treatment of Lagrange multipliers. If they are to be enforced by direct matrix transformation, then the multipoint constraints (MPCs) must have been specified using vfs_DofTabSetMPC() for each MPC equation before any element matrix information is entered into SysMatrix. SysMatrix will transform the matrix internally and eliminate all dependent degrees of freedom from the equation system. Note that transformations can not be applied to a diagonal matrix type.

To determine the structure of a sparse matrix, the SysMatrix object uses information which has been previously entered in a DofTab object. The DofTab object must be set using vfs_SysMatrixSetObject(). Once this is done use vfs_SysMatrixPreProcess() to establish the structure of the matrix. The matrix must then be initialized using vfs_SysMatrixInit(), and then assembled by calling vfs_SysMatrixAssem(), vfs_SysMatrixAssemDiag(), or vfs_SysMatrixAssemFull() for each element with the actual element matrix values. If threaded assembly is used then the assembly functions vfs_SysMatrixThreadAssem(), vfs_SysMatrixThreadAssemDiag() and vfs_SysMatrixThreadAssemFull() are used.

The system may be solved for a given right hand side vector using either a direct or iterative solution procedure. Use vfs_SysMatrixSetParami() to specify the solution procedure and preconditioner type if necessary. vfs_SysMatrixProcess() performs any pre-processing operations required before factorization occurs. vfs_SysMatrixFactor() either factorizes the matrix or computes the preconditioner depending upon the solution procedure. vfs_SysMatrixSolve() solves the system and produces a solution vector given a right hand side vector. vfs_SysMatrixSolveMult() solves the system for multiple right hand side vectors producing multiple solution vectors.

Other than boundary conditions, a mechanism is available to eliminate singularities. Use vfs_SysMatrixSetParami() with parameter SYSMATRIX_AUTOSPC set to SYS_ON. If enabled, whenever a pivot falls below the singularity threshold, the value SYSMATRIX_AUTOSPC_FACTOR multiplied by the matrix largest entry is set to the pivot to eliminate the singularity. SYSMATRIX_AUTOSPC default to SYS_OFF. SYSMATRIX_AUTOSPC_FACTOR defaults to 1.e+5.

Use vfs_SysMatrixColumn() and vfs_SysMatrixDiag() to extract a column or the diagonal of the assembled matrix into a SysVector object. Perform a matrix vector product using vfs_SysMatrixProd(). Project a matrix to a subspace defined by a set of system vectors using vfs_SysMatrixReduce().

A function can be set with vfs_SysMatrixSetFunction() to monitor the progress during a direct, sparse factorization procedure. Progress is measured by retrieving information from the SysMatrix module with vfs_SysMatrixGetInteger(), vfs_SysMatrixGetLong(), or vfs_SysMatrixGetDouble().

The factorization can be aborted with a call to vfs_SysMatrixAbort().

Destroy an instance of a SysMatrix object using vfs_SysMatrixEnd().

3.4. Linear Equation Solvers

One of the important capabilities of the SysMatrix module is its ability to solve sparse linear systems of the form Ax = b where A is an input neq by neq matrix, b is an input right hand side vector and x is an unknown vector. Two basic technologies are provided, direct sparse factorization methods and iterative preconditioned methods.

The type of solver to use depends upon the following:

  • The amount of memory available.

  • The number of equations or degrees of freedom, neq.

  • The number of right hand side vectors, b, to solve.

  • The type and polynomial order of finite elements in the model.

The direct sparse solvers are the most generally applicable to all problem sizes and types. There are two basic direct solver technologies, left looking in-core and multifrontal in-core and out-of-core. In addition there are two multifrontal solver implementations. A legacy, default method which is not BLAS optimized or parallelized and a new parallelized and BLAS optimized solver which provides for very high performance. The multi-frontal solvers may be used with in-core or out-of-core operation. A solver type is selected using vfs_SysMatrixSetParami() with type SYSMATRIX_SOLVERTYPE. To enable out-of-core storage of the factor use vfs_SysMatrixSetParami() with type SYSMATRIX_FACTOR_OOC. The file name used for the out-of-core factor may be specified using vfs_SysMatrixSetOocFile(). The out-of-core factor file is deleted when the SysMatrix object is ended using vfs_SysMatrixEnd().

Iterative PCG/GMRES solvers may provide substantial improvements in CPU time and memory. This is particularly true for 3D solid dominated models. The use of higher order elements generally improves the performance of iterative solvers. Use vfs_SysMatrixSetInitSol() to set an optional initial guess solution vector. The iterative solver preconditioner must be factorized and during iteration, back substitutions on the preconditioner are performed. Any of the linear equation solvers may be used for this purpose.

There are three convergence tolerances used in the iterative solver, relative solution error, residual error and energy error. Tolerances are specified using vfs_SysMatrixSetParamd() with parameters SYSMATRIX_ITER_UTOL, SYSMATRIX_ITER_FTOL and SYSMATRIX_ITER_ETOL respectively.

  • Solution error is assumed converged if SYSMATRIX_ITER_UTOL is achieved for two successive iterations. The solution error is evaluated in the full system.

  • Residual and energy error are assumed converged if the respective tolerances SYSMATRIX_ITER_FTOL and SYSMATRIX_ITER_ETOL are achieved in a single iteration. The residual and energy errors are evaluated in the preconditioned system.

  • Overall convergence is achieved if convergence of residual error and one of solution error or energy error is achieved.

The error measures are relative errors. If a solution or residual value used in the calculation of the relative error falls below the absolute tolerance, convergence is achieved for that error measure. This absolute tolerance is set using vfs_SysMatrixSetParamd() with parameter SYSMATRIX_ITER_ATOL. The default value is 1.E-16.

The following code fragment illustrates the basic framework of using the SysMatrix module to assemble and solve of a set of symmetric finite element matrices. The framework allows for the optional use of the iterative solver.

Vint i, j;
                /* declare objects */
vfs_SysMatrix *sysmatrix;
vfs_SysVector *sysvectorb, sysvectorx;
vfs_DofTab *doftab;
                /* number of nodes, elements */
Vint numnp, numel;
                /* number of equations, global degrees of freedom */
Vint neq;
                /* number of degrees of freedom in element,
                    element degree of freedom numbers */
Vint nedofs, lm[MAX_DOF_ELEMENT];
                /* degree of freedom types within an element */
Vint loc[MAX_DOF_ELEMENT], tag[MAX_DOF_ELEMENT];
                /* lower triangle of element stiffness by rows */
Vdouble ke[MAX_DOF_ELEMENT*(MAX_DOF_ELEMENT+1)/2];
                /* element topology and order */
Vint shape, maxi, maxj, maxk;
                /* number of nodes in element and connectivity */
Vint nix, ix[MAX_NODE_ELEMENT];
                /* node coordinates */
Vdouble x[3];

                /* create sysmatrix object */
sysmatrix = vfs_SysMatrixBegin();

doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);

for(eid = 1; eid <= numel; eid++ {
    vfs_DofTabSetMap (doftab,nedofs,loc,tag);
        /* localize element connectivity ix[] from host application */
    .
    .
    vfs_DofTabElemUse (doftab,eid,nix,ix);
}
                /* modification for iterative solver */
if(iterative_solver) {
    vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);

                /* gather and set node coordinates */
    for(j = 1; j <= numnp; j++)
        ...
        vfs_DofTabSetCoords (doftab,j,x);
    }
                /* gather and set element topology,
    for(i = 1; i <= numel; i++) {
        ...
        vfs_DofTabSetTopology (doftab,shape,maxi,maxj,maxk);
    }
    vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_ITER,neq);
} else {
    vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
}
                /* register DofTab with SysMatrix object */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab);
                /* establish matrix structure */
vfs_SysMatrixPreProcess (sysmatrix);

                /* Initialize matrix values to zero */
vfs_SysMatrixInit (sysmatrix,0.);

                /* assembly for each element */
for(i = 1; i <= numel; i++) {
                /* form element matrix and gather degrees of freedom */
    ...
    vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}

                /* equation ordering and/or symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);

                /* factorize for direct solution */
                /* form preconditioner for iterative solution */
vfs_SysMatrixFactor (sysmatrix);

                /* create load and displacement sysvector objects */
sysvectorb = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorb,neq);
sysvectorx = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorx,neq);

                /* fill load vector using
                    vfs_SysVectorAssem or
                    vfs_SysVectorScatter, etc. */

                /* solve */
vfs_SysMatrixSolve (sysmatrix,sysvectorb,sysvectorx);

                /* access displacement vector using
                    vfs_SysVectorGather */

                /* delete objects */
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (sysmatrix);
vfs_SysVectorEnd (sysvectorb);
vfs_SysVectorEnd (sysvectorx);

Note that multiple load cases may be solved by calling vfs_SysMatrixSolve() for each load case.

3.5. Eigenvalue Extraction

The SysMatrix module can also be used to solve the generalized eigenproblem. The eigenproblem may be of “standard”, buckling or vibration type. Three basic eigen solution procedures are used for large scale eigenvalue extraction, classic subspace iteration, block Lanczos and AMLS (Automatic Multi-Level Substructuring). The subspace method is used for buckling problems where a limited number of eigenvalues are to be extracted. A particular eigenvalue extraction method for vibration is selected using vfs_SysMatrixSetParami() with type SYSMATRIX_EIGEN_METHOD.

The blocked Lanczos method is specifically designed for extracting large numbers of vibration modes through a frequency range. In this case, the Lanczos method will use an automatic shifting strategy, if required, to move through the requested number of frequencies and/or frequency range from low to high frequencies. Several shifts may be required to complete the extraction of the requested vibration modes. The automatic shifting strategy is activated with any eigenvalue request which asks for the lowest modes. The initial shift in these cases is always to some point at or below the lowest frequency bound. Additional shifts are applied if all the requested eigenvalues are not found with the initial shift. Any request asking for eigenvalues nearest to a frequency will be extracted using a single shift. The user must be aware of this fact when designing the eigenvalue request.

The block size used by the Lanczos method is set using vfs_SysMatrixSetParami() with type SYSMATRIX_LANCZOS_BLOCKSIZE. In general the block size should be set to be greater than the expected eigenvalue multiplicity. For structural vibration problems with rigid body modes a value of 7 is suggested.

The AMLS eigensolution method is applicable to large scale vibration problems in which large numbers of eigenvalues are desired. The eigensolution is approximate in the sense that at successive substructures in which intermediate eigensolutions are performed, the eigensolution is truncated to some degree. The method can be several times faster than an equivalent block Lanczos solution. The AMLS method requires a factorization of the mass matrix in addition to a factorization of the stiffness matrix. Out-of-core storage of both matrices may be enabled usng vfs_SysMatrixSetParami() with type SYSMATRIX_FACTOR_OOC. The names of the out-of-core files may be specified using vfs_SysMatrixSetOocFile(). The degree of eigenvalue truncation affects the accuracy and runtime of the AMLS method. Use SYSMATRIX_AMLS_CUTFACT1 to set the frequency cutoff factor. This factor is multiplied by the high frequency limit of the requested frequency spectrum. Use SYSMATRIX_AMLS_CUTFACT2 to control a secondary cutoff frequency used in the reduced problem. This parameter is particularly effective in increasing the accuracy of the highest frequencies.

The subspace method should be used only for buckling problems. There is no automatic shifting strategy available for this method due to the normally low number of modes requested.

The function vfs_SysMatrixEigen() is designed to compute a subset of eigenvalues and eigenvectors for large systems using the subspace, Lanczos or AMLS method. The function vfs_SysMatrixEigenAll() is designed to compute all eigenvalues and eigenvectors for the generalized eigenproblem for small systems.

The computation of eigenvectors is performed during the eigenvalue extraction process for the subspace and block Lanczos eigensolvers. It is performed, on demand, as a post processing step for the AMLS eigensolver. Therefore, an option exists for out-of-core storage of the eigenvectors to reduce memory usage for the subspace and Lanczos eigensolvers, use vfs_SysMatrixSetParami() with type SYSMATRIX_EIGEN_OOC. The file name used for the out-of-core vectors may be specified using vfs_SysMatrixSetOocFile(). Use the function vfs_SysMatrixEigenMode() to compute/retrieve eigenvectors. For the subspace and Lanczos methods, the pre-computed eigenvectors are retrieved from memory or read from the specified out-of-core file. For the AMLS method the eigenvectors are computed when requested in memory blocks for computational efficiency. The requested eigenvector is selected from the currently computed block. If there are memory constraints, use the function vfs_SysMatrixSetParami() with type SYSMATRIX_AMLS_RECOVERSIZE to set the number of vectors in the computational block. Since eigenvectors are computed on-demand for the AMLS method, the vfs_SysMatrixEigenMode() function call can take considerable time if the computation of a new block of eigenvectors is triggered by the call. After all desired eigenvectors are retrieved, free memory and delete any out-of-core file using vfs_SysMatrixEigenFree().

3.6. Function Descriptions

The currently available SysMatrix functions are described in detail in this section.

vfs_SysMatrix *vfs_SysMatrixBegin(void)

create an instance of a SysMatrix object

Create an instance of a SysMatrix object. Memory is allocated for the object private data and the pointer to the object is returned.

Destroy an instance of a SysMatrix object using 

void vfs_SysMatrixEnd (vfs_SysMatrix *sysmatrix)

Return the current value of a SysMatrix object error flag using 

Vint vfs_SysMatrixError (vfs_SysMatrix *sysmatrix)

Returns:

The function returns a pointer to the newly created SysMatrix object. If the object creation fails, NULL is returned.

void vfs_SysMatrixEnd(vfs_SysMatrix *p)

destroy an instance of a SysMatrix object

See vfs_SysMatrixBegin()

Vint vfs_SysMatrixError(vfs_SysMatrix *p)

return the current value of a SysMatrix object error flag

See vfs_SysMatrixBegin()

void vfs_SysMatrixDef(vfs_SysMatrix *p, Vint type, Vint neq)

define number of global dof and structure

Specify the number of global dof and structure of the matrix. Any previously assembled or factored system matrices are deleted.

Inquire of a defined type and neq as output arguments using 

void vfe_SysMatrixInq (vfe_SysMatrix *sysmatrix,
                       Vint *type,
                       Vint *neq)

Errors

  • SYS_ERROR_VALUE is generated if neq is less than or equal to zero.

  • SYS_ERROR_ENUM is generated if an improper type is input.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Type of matrix 

    x=SYSMATRIX_SYMM_SPARSE     Symmetric sparse
 =SYSMATRIX_NSYMM_SPARSE    Non-symmetric sparse
 =SYSMATRIX_SYMM_FULL       Symmetric full
 =SYSMATRIX_NSYMM_FULL      Non-symmetric full
 =SYSMATRIX_SYMM_ITER       Symmetric iterative
 =SYSMATRIX_NSYMM_ITER      Non-symmetric iterative
 =SYSMATRIX_DIAG            Diagonal

  • neq – Number of global dof

void vfs_SysMatrixInit(vfs_SysMatrix *p, Vdouble v)

initialize matrix to a given value

Initialize all values in the assembled system matrix to v.

Parameters:

  • p – Pointer to SysMatrix object.

  • v – Initialization value

void vfs_SysMatrixInq(vfs_SysMatrix *p, Vint *type, Vint *neq)

define number of global dof and structure

See vfs_SysMatrixDef()

void vfs_SysMatrixSetParami(vfs_SysMatrix *p, Vint type, Vint param)

set solution parameters

Set solution parameters.

Specify the eigenproblem procedure using SYSMATRIX_EIGEN. By default SYSMATRIX_EIGEN is SYSMATRIX_EIGEN_STAN.

Specify the number of eigenvalues sought using SYSMATRIX_EIGEN_NUM. The number of requested eigenvalues must be less than or equal to the number of global dof set using vfs_SysMatrixDef(). The default is one.

Set the eigensolution type using SYSMATRIX_EIGEN_TYPE. This parameter determines whether the nearest eigenvalues about a given shift are sought, SYS_EIGEN_NEAREST; the lowest eigenvalues in an interval are sought, SYS_EIGEN_LOWEST; or all the eigenvalues in an interval are sought, SYS_EIGEN_ALL. In any case the eigensolution will terminate if the specified number of eigenvalues set using the SYSMATRIX_EIGEN_NUM parameter are computed. The interval limits are set using vfs_SysMatrixSetParamd(). By default SYSMATRIX_EIGEN_TYPE is SYS_EIGEN_LOWEST.

Specify the maximum number of eigenvalue iterations taken at each shift using SYSMATRIX_EIGEN_MAXITER. The default is 50.

Specify the normalization of eigenvectors using SYSMATRIX_EIGEN_NORM. If the parameter value is SYS_EIGEN_NORMMAX, then the eigenvector is normalized to unit value of the largest component. If the parameter value is SYS_EIGEN_NORMMASS, then the eigenvector is normalized to unit value of the generalized mass. By default SYSMATRIX_EIGEN_NORM is SYS_EIGEN_NORMMAX.

Specify the existence of rigid body modes (zero eigenvalues) using SYSMATRIX_EIGEN_ZEROMODES. If it is known that no rigid body modes exist in the eigensolution then it is computationally advantageous to disable this flag. By default SYSMATRIX_EIGEN_ZEROMODES is SYS_ON.

Use SYSMATRIX_EIGEN_OOC to indicate that eigenvectors are to be stored out of core.

The parameter SYSMATRIX_FREEBLASBUFFERS toggles the freeing of external BLAS library buffers. By default SYSMATRIX_FREEBLASBUFFERS is set to SYS_OFF.

The toggle SYSMATRIX_SKIPZERODIAG allows for the factorization and solution of a diagonal matrix where some entries are zero. If set to SYS_ON then zero entries will be skipped and the matrix is assumed non-singular. Defaults to SYS_OFF.

The parameter SYSMATRIX_CACHESIZE sets the central processor cache size in kilobytes. This value should be set to the hardware cache size as closely as possible for optimum performance of factor and solve operations. By default SYSMATRIX_CACHESIZE is 256.

The parameter SYSMATRIX_WORKSIZE_FACTOR sets the block size in megabytes to be used with the out of core solver option. By default SYSMATRIX_WORKSIZE_FACTOR is 256.

Use SYSMATRIX_FACTOR_OOC to indicate that an out-of-core solver is to be used.

The toggle SYSMATRIX_AUTOSPC enables an algorithm that penalizes a detected singularity so that factorization can proceed. This parameter should be used in conjunction with the floating point parameters SYSMATRIX_FACTOR_SINGTOL and SYSMATRIX_AUTOSPC_FACTOR. By default SYSMATRIX_AUTOSPC is set to SYS_OFF.

Use SYSMATRIX_MONITOR_FREQ to specify the frequency to be used to call a monitoring function during a sparse factorization procedure. Defaults to 1000.

Use SYSMATRIX_EIGEN_INITVECTYPE to set the algorithm that will be employed during the initialization phase of a vibration problem. If set to 0, the first eigenvector is set equal to the scaled diagonal of the mass matrix, while all subsequent vectors are randomly initialized with values between -1.0 and 1.0. If set to 1, then the first vector a random value between -0.1 and 0.1 is added to all degrees of freedom for which the diagonal of the mass matrix is larger than a tolerance specified with vfs_SysMatrixSetParamd() with the parameter SYSMATRIX_EIGEN_MASSTOL. An array is then constructed with the ratio of the diagonal of the mass and stiffness matrices is computed. If the stiffness matrix entry is smaller than the parameter SYSMATRIX_EIGEN_STIFFTOL then the ratio for the degree of freedom is set to zero. This array is sorted. For the second vector, the degree of freedom associated with the largest ratio is set to 1.0, while all others for which the mass diagonal is larger than SYSMATRIX_EIGEN_MASSTOL are set to a random value between -0.1 and 0.1. For the third vector, the degree of freedom with the second largest ratio is chosen and the algorithm is repeated for all other vectors. Defaults to 0.

Use SYSMATRIX_EIGEN_METHOD to set the eigensolver type used for vibration analysis in vfs_SysMatrixEigen(). Method can be SYS_EIGEN_SUBSPACE, SYS_EIGEN_LANCZOS or SYS_EIGEN_AMLS. Defaults to SYS_EIGEN_LANCZOS. The SYS_EIGEN_AMLS method can only be used if the solver type is SYS_SOLVERTYPE_MFP. Note the block Lanczos and AMLS eigensolvers require the MKL BLAS or OpenBLAS.

Use SYSMATRIX_AMLS_RECOVERSIZE to set the eigenvector block size used during the eigenvector recovery computation in the AMLS eigensolver. This should only be changed to a smaller number if there are memory constraints. The default is set to 128.

Use SYSMATRIX_LANCZOS_BLOCKSIZE to select the block size used by the Lanczos eigensolver. If set to 0 a simplified unblocked Lanczos algorithm is used which avoids the need for an external BLAS library. Note that performance is significantly degraded if an external BLAS library is not used. By default, if VKI_LIBAPI_BLASMKL_SEQUENTIAL, VKI_LIBAPI_BLASMKL_THREAD or VKI_LIBAPI_OPENBLAS is defined, then SYSMATRIX_LANCZOS_BLOCKSIZE defaults to 7, else it defaults to 0. An error will be issued if the block size is set to a value greater than zero and some version of the MKL BLAS has not been defined. Generally the block size should be greater than or equal to the largest expected eigenvalue multiplicity. For structural vibration problems with rigid body modes a value of 7 is appropriate.

Use SYSMATRIX_TECH to set the linear equation technology to be used. The technology defaults to SYSMATRIX_TECH_SPARSE for any of the direct methods, to SYSMATRIX_TECH_CONJGRAD for SYSMATRIX_SYMM_ITER, and to SYSMATRIX_TECH_GMRES for SYSMATRIX_NSYMM_ITER. Use SYSMATRIX_TECH_GMRES for matrices of type SYSMATRIX_SYMM_ITER if the matrix is expected to be indefinite.

Use SYSMATRIX_SOLVERTYPE to set the specific direct sparse equation solver package to be used. Set SYS_SOLVERTYPE_LL to use the serial left-looking, in-core only solver, SYS_SOLVERTYPE_MF to use the legacy serial multi-frontal solver and SYS_SOLVERTYPE_MFP to use the parallel multi-frontal solver. Note the parallel multi-frontal solver requires the MKL BLAS. Set SYS_SOLVERTYPE_PARDISO to use Intel’s MKL Pardiso solver. Pardiso selection requires that DevTools be compiled with VKI_LIBAPI_PARDISO. Set SYS_SOLVERTYPE_MUMPS to use the NUMPS solver. MUMPS selection requires that DevTools be compiled with VKI_LIBAPI_MUMPS. Defaults to SYS_SOLVERTYPE_MFP.

Errors

  • SYS_ERROR_ENUM is generated if an improper type is specified.

  • SYS_ERROR_ENUM is generated if an improper param is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Integer parameter type to set. 

    x=SYSMATRIX_EIGEN              Eigenproblem procedure
 =SYSMATRIX_EIGEN_LOWERTEST    Test for lower bound on eigenvalue
 =SYSMATRIX_EIGEN_UPPERTEST    Test for upper bound on eigenvalue
 =SYSMATRIX_EIGEN_MAXITER      Maximum iterations per shift
 =SYSMATRIX_EIGEN_METHOD       Eigensolver method for vibration
 =SYSMATRIX_EIGEN_NORM         Eigenvector normalization
 =SYSMATRIX_EIGEN_NUM          Number of eigenvalues requested
 =SYSMATRIX_EIGEN_TYPE         Eigensolution type
 =SYSMATRIX_EIGEN_ZEROMODES    Rigid body modes flag
 =SYSMATRIX_EIGEN_INITVECTYPE  Vibration eigenvector initialization type
 =SYSMATRIX_EIGEN_OOC          Out-of-core eigenvector flag
 =SYSMATRIX_FREEBLASBUFFERS    Free external BLAS library buffers
 =SYSMATRIX_AMLS_RECOVERSIZE   AMLS vector recovery block size
 =SYSMATRIX_LANCZOS_BLOCKSIZE  Lanczos block size
 =SYSMATRIX_TECH               Solution procedure
 =SYSMATRIX_CACHESIZE          Cache size (KiloBytes)
 =SYSMATRIX_WORKSIZE_FACTOR    Out of core factor memory (MB)
 =SYSMATRIX_SOLVETEMP          Toggle applied displacements during solve
 =SYSMATRIX_MAXCGITER          Max iterative solver iterations
 =SYSMATRIX_SKIPZERODIAG       Toggle to ignore zero entries in diagonal matrix
 =SYSMATRIX_FACTOR_OOC         Out-of-core matrix factorization flag
 =SYSMATRIX_AUTOSPC            Toggle to penalize singularities
 =SYSMATRIX_MONITOR_FREQ       Monitoring frequency for factorization
 =SYSMATRIX_SOLVERTYPE         Linear equation solver

  • param – Specifies the integer value that type will be set to. 

    x=SYSMATRIX_EIGEN_BUCK         Buckling eigenproblem
 =SYSMATRIX_EIGEN_GENERAL      Generalized eigenproblem
 =SYSMATRIX_EIGEN_RIGID        Rigid body eigenproblem
 =SYSMATRIX_EIGEN_STAN         Standard eigenproblem
 =SYSMATRIX_EIGEN_VIBE         Vibration eigenproblem
 =SYS_EIGEN_AMLS               AMLS method for vibration
 =SYS_EIGEN_LANCZOS            Lanczos method for vibration
 =SYS_EIGEN_SUBSPACE           Subspace method for vibration
 =SYS_EIGEN_ALL                Get all eigenvalues in range
 =SYS_EIGEN_LOWEST             Get lowest eigenvalues in range
 =SYS_EIGEN_NEAREST            Get nearest eigenvalues
 =SYS_EIGEN_NORMMASS           Mass normalization
 =SYS_EIGEN_NORMMAX            Max component normalization
 =SYS_ON                       Enable
 =SYS_OFF                      Disable
 =SYS_SOLVERTYPE_LL            Left looking solver
 =SYS_SOLVERTYPE_MF            Multifrontal solver
 =SYS_SOLVERTYPE_MFP           Multifrontal parallel solver
 =SYS_SOLVERTYPE_PARDISO       Intel Pardiso solver
 =SYS_SOLVERTYPE_MUMPS         MUMPS solver

void vfs_SysMatrixSetParamd(vfs_SysMatrix *p, Vint type, Vdouble param)

set solution parameters

Set solution parameters. If the iterative solver solution procedure is specified then the relative convergence tolerance for displacements may be set with SYSMATRIX_ITER_UTOL By default SYSMATRIX_ITER_UTOL is set to 5.E-7. The displacement convergence measure is the maximum displacement change component over the total displacement at the same component. The relative convergence tolerance for force is specified with SYSMATRIX_ITER_FTOL. The force convergence measure is the L2 norm of the residual over the L2 norm of the force. By default SYSMATRIX_ITER_FTOL is set to 5.E-3. The relative convergence tolerance for energy is specified with SYSMATRIX_ITER_ETOL. By default SYSMATRIX_ITER_ETOL is set to 5.E-7. The absolute convergence tolerance for displacement and force is specified with SYSMATRIX_ITER_ATOL. By default SYSMATRIX_ITER_ATOL is set to 1.E-16.

Specify the lower limit of the eigenvalue interval using SYSMATRIX_EIGEN_LOWER and the upper limit using SYSMATRIX_EIGEN_UPPER. The defaults are [0.,1.].

Specify a given eigenvalue shift using SYSMATRIX_EIGEN_SHIFT. The default is 0.

Specify an overall cutoff factor at substructure levels using SYSMATRIX_AMLS_CUTFACT1. The frequencies at each substructure level are truncated at this factor times the upper frequency limit of requested spectrum. This value will affect the accuracy of all modes. Suggest specified values be no larger than 10. By default SYSMATRIX_AMLS_CUTFACT1 is 5.

Specify a reduced cutoff factor using SYSMATRIX_AMLS_CUTFACT2. Higher values will have a tendency to increase the accuracy of higher modes. By default SYSMATRIX_AMLS_CUTFACT2 is 1.7 Suggest specified values be no larger than 2.5. Note that 1. < SYSMATRIX_AMLS_CUTFACT2 < SYSMATRIX_AMLS_CUTFACT1.

Specify a subspace eigenvalue convergence tolerance using SYSMATRIX_SUBSPACE_UTOL. The default is .0001 . Specify a Lanczos eigenvalue convergence tolerance using SYSMATRIX_LANCZOS_UTOL. The default is 3.7e-11 It is not recommended to alter these values.

Specify a factorization singularity tolerance using SYSMATRIX_FACTOR_SINGTOL. This tolerance is used to compare the original diagonal pivot to its reduced value. If the integer parameter SYSMATRIX_AUTOSPC is set to SYS_ON, if a singularity is encountered then the reduced pivot is replaced by the value given by the parameter SYSMATRIX_AUTOSPC_FACTOR multiplied by the matrix largest entry. This procedure is equivalent to introducing a penalty parameter that multiplies the matrix diagonal in order to enforce a constraint. Users should exert care in using this procedure, especially if there is a load applied in the direction of the singularity, or if the system is indefinite. By default SYSMATRIX_FACTOR_SINGTOL is 1.E-13 and SYSMATRIX_AUTOSPC_FACTOR is 1.e+5.

Specify an absolute singularity tolerance using SYSMATRIX_FACTOR_ABSZERO. The system is considered singular if a pivot’s absolute value falls below this threshold.

Errors

  • SYS_ERROR_ENUM is generated if an improper type is specified.

  • SYS_ERROR_VALUE is generated if an improper param is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Parameter type to set. 

    x=SYSMATRIX_ITER_UTOL       Convergence tolerance displacement
 =SYSMATRIX_ITER_FTOL       Convergence tolerance force
 =SYSMATRIX_ITER_ETOL       Convergence tolerance energy
 =SYSMATRIX_ITER_ATOL       Convergence tolerance absolute
 =SYSMATRIX_EIGEN_LOWER     Lower limit of eigenvalue interval
 =SYSMATRIX_EIGEN_UPPER     Upper limit of eigenvalue interval
 =SYSMATRIX_EIGEN_SHIFT     Given eigenvalue shift
 =SYSMATRIX_AMLS_CUTFACT1   Overall frequency cutoff factor
 =SYSMATRIX_AMLS_CUTFACT2   Reduced frequency cutoff factor
 =SYSMATRIX_LANCZOS_UTOL    Lanczos eigenvalue tolerance
 =SYSMATRIX_SUBSPACE_UTOL   Subspace eigenvalue tolerance
 =SYSMATRIX_AUTOSPC_FACTOR  Factor singularity penalty multiplier
 =SYSMATRIX_FACTOR_SINGTOL  Relative singularity tolerance
 =SYSMATRIX_FACTOR_ABSZERO  Absolute singularity tolerance

  • param – Specifies the double value that type will be set to.

void vfs_SysMatrixSetOocFile(vfs_SysMatrix *p, Vint type, Vchar *file)

set out-of-core filename

Set out-of-core file name. These files are temporary files that are created and deleted by the SysMatrix object. Use SYSMATRIX_KFACTOR_FILE to specify the name of the file to be used to store the stiffness matrix factor during an out-of-core matrix factorization. Defaults to vki.foc.

Use SYSMATRIX_MFACTOR_FILE to specify the name of the file to be used to store the mass matrix factor during an out-of-core matrix factorization. Defaults to vki.moc. A mass matrix factor is generated by the AMLS eigensolution procedure.

Use SYSMATRIX_EIGEN_FILE to specify the name of the file to be used to store the eigenvectors during eigen extraction. Defaults to vki.eoc.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Parameter type to set. 

    x=SYSMATRIX_KFACTOR_FILE    Out-of-core stiffness factor
 =SYSMATRIX_MFACTOR_FILE    Out-of-core mass factor
 =SYSMATRIX_EIGEN_FILE      Out-of-core eigenvectors

  • file – Specifies the name of the file to be used

void vfs_SysMatrixSetValue(vfs_SysMatrix *p, Vint i, Vint j, Vdouble value[])

set matrix value

Set the value of an individual entry in the matrix. For real matrices, only the first entry of value is used. For complex matrices, both values are used. vfs_SysMatrixSetComplexMode() or vfs_SysMatrixSetThreadComplexMode() do not apply; both real and imaginary parts must be defined for complex matrices.

The column and row values, i and j, must come from the matrix structure specified by the user with vfs_DofTabPointersIndicesPtr(). SysMatrix will ensure that, if any reordering is performed, the specified value will be assigned to the properly permuted matrix.

Errors

SYS_ERROR_OPERATION is generated if entry a_ij is not found or if vfs_DofTabPointersIndicesPtr() was not called.

Parameters:

  • p – Pointer to SysMatrix object.

  • i – 0-based column number.

  • j – 0-based row number.

  • value – Value to set matrix entry to.

void vfs_SysMatrixAssem(vfs_SysMatrix *p, Vint nedofs, Vint lm[], Vdouble s[])

assemble an element matrix

Use vfs_SysMatrixAssem() to assemble the element matrix s in lower triangular form, The number of element degrees of freedom for the element is nedofs. The global dof numbers are input in the lm vector.

Errors

  • SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

  • SYS_ERROR_OPERATION is generated if the matrix has not been initialized with vfs_SysMatrixInit().

Parameters:

  • p – Pointer to SysMatrix object.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • s – Element matrix

void vfs_SysMatrixThreadAssem(vfs_SysMatrix *p, Vint ithread, Vint nedofs, Vint lm[], Vdouble s[])

assemble an element matrix

Assemble an element matrix on thread ithread. These functions are otherwise identical to their serial counterparts, vfs_SysMatrixAssem(), vfs_SysMatrixAssemDiag() and vfs_SysMatrixAssemFull().

Errors

  • SYS_ERROR_OPERATION is generated if vfs_SysMatrixSetNumThreadsAssem() has not been previously called.

  • SYS_ERROR_VALUE is generated if ithread is less than 1 or greater than the maximum number of threads specified in vfs_SysMatrixSetNumThreadsAssem().

  • SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

  • SYS_ERROR_OPERATION is generated if the matrix has not been initialized with vfs_SysMatrixInit().

Parameters:

  • p – Pointer to SysMatrix object.

  • ithread – Thread number starting from 1.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • s – Element matrix

void vfs_SysMatrixAssemFull(vfs_SysMatrix *p, Vint nedofs, Vint lm[], Vdouble sf[])

assemble a full element matrix

Use vfs_SysMatrixAssemFull() to assemble a full matrix. The number of element degrees of freedom for the element is nedofs. The global dof numbers are input in the lm vector.

Errors

  • SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

  • SYS_ERROR_OPERATION is generated if the matrix has not been initialized with vfs_SysMatrixInit().

Parameters:

  • p – Pointer to SysMatrix object.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • sf – Element full matrix

void vfs_SysMatrixAssemDiag(vfs_SysMatrix *p, Vint nedofs, Vint lm[], Vdouble sd[])

assemble a element matrix in diagonal form

Use vfs_SysMatrixAssemDiag() to assemble the element matrix sd in diagonal form. The number of element degrees of freedom for the element is nedofs. The global dof numbers are input in the lm vector.

Errors

  • SYS_ERROR_VALUE is generated if an improper nedofs or lm entry is specified.

  • SYS_ERROR_OPERATION is generated if the matrix has not been initialized with vfs_SysMatrixInit().

Parameters:

  • p – Pointer to SysMatrix object.

  • nedofs – Number of element degrees of freedom

  • lm – Vector of element global dof numbers

  • sd – Element diagonal matrix

void vfs_SysMatrixFactor(vfs_SysMatrix *p)

factorize

Factorize a system matrix. If a direct solver has been specified, a complete factorization is performed. If an iterative solver has been specified, then the preconditioner is computed. The function vfs_SysMatrixProcess() must have been previously called.

Errors

SYS_ERROR_OPERATION is generated if the matrix has not been assembled or has not been processed.

Parameters:

p – Pointer to SysMatrix object.

void vfs_SysMatrixSolve(vfs_SysMatrix *p, vfs_SysVector *sysvectorx, vfs_SysVector *sysvectory)

solve given a system vector

Solve for system vector sysvectory given system vector sysvectorx. If the flag SYSMATRIX_SOLVETEMP is set, values in sysvectory for which there is a temporary restraint set with vfs_SysMatrixSetTemp() are taken as applied boundary conditions; the given right-hand-side vector sysvectorx will be modified to account for these values.

Errors

  • SYS_ERROR_OPERATION is generated if the number of global dof in each system matrix and system vector are not equal.

  • SYS_ERROR_COMPUTE is generated if, during an iterative solution, either it is detected that the matrix is not positive definite, or the maximum number of iterations is reached.

Parameters:
void vfs_SysMatrixSolveMult(vfs_SysMatrix *p, Vint nvec, vfs_SysVector *sysvectorx[], vfs_SysVector *sysvectory[])

solve given an array of system vectors

Solve for nvec system vectors sysvectory given system vectors sysvectorx. This is a multiple vector version of vfs_SysMatrixSolve().

Errors

  • SYS_ERROR_OPERATION is generated if the number of global dof in each system matrix and system vector are not equal.

  • SYS_ERROR_COMPUTE is generated if, during an iterative solution, either it is detected that the matrix is not positive definite, or the maximum number of iterations is reached

Parameters:

  • p – Pointer to SysMatrix object.

  • nvec – Number of system vectors 0 < nvec

  • sysvectorx – Pointer to array of nvec SysVector objects

  • sysvectory – Pointer to array of nvec SysVector objects

void vfs_SysMatrixSetTemp(vfs_SysMatrix *p, Vint ndofs, Vint lm[])

set temporary restraints

Set ndofs global dof numbers lm as temporary restraints.

Errors

SYS_ERROR_VALUE is generated if an improper ndofs is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • ndofs – Number of degrees of freedom

  • lm – Vector of global dof numbers to restrain

void vfs_SysMatrixColumn(vfs_SysMatrix *p, Vint lmcol, vfs_SysVector *sysvector)

extract a column

Extract a column of sysmatrix and place it in sysvector. The column number is in the interval 1 <= lmcol <= number of global dof. Use vfs_SysMatrixDiag() to extract the diagonal.

Errors

  • SYS_ERROR_OPERATION is generated if the number of global dof in the system matrix and system vector is not equal.

  • SYS_ERROR_VALUE is generated if an improper lmcol is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • lmcol – Column number to extract.

  • sysvector – Pointer to SysVector object

void vfs_SysMatrixProd(vfs_SysMatrix *p, vfs_SysVector *sysvectorx, vfs_SysVector *sysvectory)

multiply by a system vector

Compute the product of system matrix sysmatrix with system vector sysvectorx and place the result in system vector sysvectory.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system matrix and system vector are not equal.

Parameters:
void vfs_SysMatrixDiag(vfs_SysMatrix *p, vfs_SysVector *sysvector)

extract the diagonal

Extract the diagonal of sysmatrix and place it in sysvector. Use vfs_SysMatrixColumn() to extract a column.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in the system matrix and system vector is not equal.

Parameters:
void vfs_SysMatrixReduce(vfs_SysMatrix *p, Vint nvec, Vint nattach, Vint nvibe, Vint lmattach[], vfs_SysVector *sysvector[], vfs_SysMatrix *sysmatrixq)

reduce a matrix

Reduce the system matrix, sysmatrix by the nvec system vectors, sysvector[], and place the result in the system matrix sysmatrixq. The reduced system matrix, sysmatrixq must be a full matrix of order equal to nvec.

If the matrix to be reduced is the stiffness matrix then there is a considerable computational efficiency in specifying those modes which are attachment modes. Attachment modes are those modes due to applied unit displacement at a specific degree of freedom. These are the “retained” degrees of freedom. These nattach modes must appear first in the set of sysvector modes and nattach related global degrees of freedom must be specified in lmattach. If vibration modes are present they must come after the attachment points. If the matrix to be reduced is not a stiffness matrix (eg. mass matrix), set nattach to zero and lmattach may be input NULL.

In matrix form, given 

sysvector[0] = [ a_00 a_01 ... a_0N ]
sysvector[1] = [ a_10 a_11 ... a_1N ]
...
sysvector[n] = [ a_n0 a_n1 ... a_nN ]
then 
             --                  --           --                  --
             | a_00 a_01 ... a_0N |           | a_00 a_10 ... a_n0 |
sysmatrixq = | a_10 a_11 ... a_1N | sysmatrix | a_01 a_11 ... a_n1 |
             |     ...            |           |         ...        |
             | a_n0 a_n1 ... a_nN |           | a_0N a_1N ... a_nN |
             --                  --           __                  __
Errors

  • SYS_ERROR_VALUE is generated if an improper nvec is specified.

  • SYS_ERROR_OPERATION is generated if the number of global dof in sysmatrix and sysvector[] are not equal.

  • SYS_ERROR_OPERATION is generated if sysmatrixq is not a full matrix.

Parameters:

  • p – Pointer to SysMatrix object.

  • nvec – Number of system vectors 0 < nvec

  • nattach – Number of attachment modes.

  • nvibe – Number of vibration modes.

  • lmattach – Vector of attachment mode global dof numbers

  • sysvector – Pointer to array of nvec SysVector objects

  • sysmatrixq – Pointer to SysMatrix object of reduced matrix

void vfs_SysMatrixInner(vfs_SysMatrix *p, Vint nvec, vfs_SysVector *sysvectorx[], vfs_SysVector *sysvectory[])

form matrix as inner product of vectors

Form the matrix A as the inner product, A=X(transposed)*Y, of the set of nvec vectors X and Y. The matrix A must be a symmetric full matrix of order nvec.

Errors

  • SYS_ERROR_VALUE is generated if the matrix A is not of order nvec.

  • SYS_ERROR_OPERATION is generated if the matrix A is not a full symmetric matrix.

Parameters:

  • p – Pointer to SysMatrix object.

  • nvec – number of vectors X and Y.

  • sysvectorx – Pointer to vector of SysVector objects X

  • sysvectory – Pointer to vector of SysVector objects Y

void vfs_SysMatrixEigenAll(vfs_SysMatrix *p, vfs_SysMatrix *sysmatrixb, vfs_SysVector *sysvectore, vfs_SysVector *sysvectorx[])

extract all eigenvalues and eigenvectors

Solve the generalized eigenproblem (A-lambda*B)X = 0 . The system matrices must be in SYSMATRIX_SYMM_FULL format. Since all the eigenvalues and eigenvectors are to be solved for, this procedure is only practical when the number of global dof is on the order of 100’s.

The type of eigenproblem solved is dependent upon the SYSMATRIX_EIGEN parameter set using vfs_SysMatrixSetParami(). If the eigenproblem type is SYSMATRIX_EIGEN_STAN then the matrix B is assumed to be the identity and sysmatrixb is ignored. If the eigenproblem type is SYSMATRIX_EIGEN_GENERAL then the matrix B must be input and the eigenvalues as defined above are returned. If the eigenproblem type is SYSMATRIX_EIGEN_BUCK then the matrix B is assumed to be the geometric stiffness matrix. The eigenvalues returned in sysvectore will be set to - lambda. If the eigenproblem type is SYSMATRIX_EIGEN_VIBE then the matrix B is assumed to be a mass matrix. The eigenvalues returned in sysvectore will be set to sqrt(lambda)/(2*pi). If the eigenproblem type is SYSMATRIX_EIGEN_RIGID then the matrix B is assumed to be the identity and sysmatrixb is ignored.

Errors

SYS_ERROR_OPERATION is generated if sysmatrix or sysmatrixb is not defined with a SYSMATRIX_SYMM_FULL structure.

Parameters:

  • p – Pointer to SysMatrix object.

  • sysmatrixb – Pointer to SysMatrix object for matrix B.

  • sysvectore – Pointer to SysVector object of eigenvalues

  • sysvectorx – Pointer to array of SysVector objects of eigenvectors X.

void vfs_SysMatrixEigen(vfs_SysMatrix *p, vfs_SysMatrix *sysmatrixb, Vint *neigen, Vdouble eigenval[])

extract set of eigenvalues and eigenvectors

Solve the generalized eigenproblem (A-lambda*B)X = 0 for a set of eigenvalues. The number of eigenvalues computed is returned in neigen and the computed eigenvalues are returned in eigenval. Retrieve the eigenvectors using vfs_SysMatrixEigenMode(). Free any storage associated with the internal storage of eigenvectors using vfs_SysMatrixEigenFree().

The eigenvalue/eigenvector pairs computed by this function depend on the values of integer and double precision parameters set, respectively, with vfs_SysMatrixSetParami() and vfs_SysMatrixSetParamd(). The supported combinations are outlined in the table below.

Supported eigenvalue parameter combinations

Case

SYSMATRIX_EIGEN_TYPE

SYSMATRIX_EIGEN_LOWERTEST

SYSMATRIX_EIGEN_UPPERTEST

SYSMATRIX_EIGEN_NUM

1

SYS_EIGEN_ALL

SYS_ON

SYS_ON

N/A

2

SYS_EIGEN_ALL

SYS_OFF

SYS_ON

N/A

3

SYS_EIGEN_LOWEST

SYS_ON

SYS_ON

0

4

SYS_EIGEN_LOWEST

SYS_ON

SYS_ON

n

5

SYS_EIGEN_LOWEST

SYS_OFF

SYS_ON

0

6

SYS_EIGEN_LOWEST

SYS_OFF

SYS_ON

n

7

SYS_EIGEN_LOWEST

SYS_ON

SYS_OFF

0

8

SYS_EIGEN_LOWEST

SYS_ON

SYS_OFF

n

9

SYS_EIGEN_NEAREST

N/A

N/A

0

10

SYS_EIGEN_NEAREST

N/A

N/A

n

Case 1 will compute all eigenvalues found in the interval defined by [SYSMATRIX_EIGEN_LOWER,SYSMATRIX_EIGEN_UPPER]. The number of requested eigenvalues is ignored, and the user should ensure that the eigenval argument is sufficiently large to hold all eigenvalues encountered in the interval.

Case 2 is similar to case 1, except that the interval is replaced by [0,SYSMATRIX_EIGEN_UPPER].

Case 3 computes the single lowest eigenvalue in the interval [SYSMATRIX_EIGEN_LOWER,SYSMATRIX_EIGEN_UPPER].

Case 4 is similar to Case 3, except that the n lowest eigenvalues in the interval are computed.

Cases 5 and 6 are similar to Cases 3 and 4, respectively, except that the interval is replaced by [0,SYSMATRIX_EIGEN_UPPER].

Cases 7 and 8 are similar to cases 3 and 4, respectively, except that the interval is replaced by [SYSMATRIX_EIGEN_LOWER,infinity].

Case 9 computes the single eigenvalue closest to the shift SYSMATRIX_EIGEN_SHIFT, while Case 10 computes the n closest eigenvalues to the shift SYSMATRIX_EIGEN_SHIFT.

When using the AMLS eigensolver, only Case 4 and Case 6 are supported. The number of eigenvalues and upper bound must be specified. The lower bound is optional.

Errors

SYS_ERROR_COMPUTE is generated if sysmatrixb has all zero diagonals.

Parameters:

  • p – Pointer to SysMatrix object.

  • sysmatrixb – Pointer to SysMatrix object for matrix B.

  • neigen – [out] Number of eigenvalues found

  • eigenval – [out] Vector of neigen computed eigenvalues

void vfs_SysMatrixEigenMode(vfs_SysMatrix *p, Vint ieigen, vfs_SysVector *sysvector)

retrieve a computed eigenvector

Return the eigenvector of computed eigenmode ieigen in sysvector. The eigen mode number must be in the range 1 <= ieigen <= neigen where neigen is the number of eigenvalues returned by vfs_SysMatrixEigen().

Errors

SYS_ERROR_VALUE is generated if an improper ieigen is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • ieigen – Eigen mode to retrieve

  • sysvector – Pointer to SysVector for eigenvector

void vfs_SysMatrixEigenFree(vfs_SysMatrix *p)

free all computed eigenvector storage

Free all storage used to maintain eigenvectors computed by vfs_SysMatrixEigen(). Call this function after all desired eigenvectors have been retrieved using vfs_SysMatrixEigenMode().

Parameters:

p – Pointer to SysMatrix object.

void vfs_SysMatrixSum(vfs_SysMatrix *p, Vdouble s, vfs_SysMatrix *sysmatrixb)

sum a system matrix

Sum s times sysmatrixb to sysmatrix. The operation replaces any previous values in sysmatrix with the sum of sysmatrix and s times sysmatrixb.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in each system matrix are not equal or if the matrix types are dissimilar.

Parameters:

  • p – Pointer to SysMatrix object.

  • s – number of sums

  • sysmatrixb – Pointer to SysMatrix object b

void vfs_SysMatrixGetInteger(vfs_SysMatrix *p, Vint type, Vint param[])

get integer matrix information

Query for integer information associated with the system matrix.

The integer queries SYSMATRIX_NUMNEG and SYSMATRIX_EXPONENT require the system matrix to have been factorized using vfs_SysMatrixFactor(). SYSMATRIX_MAXSUPERSIZE is the size of the largest supernode used during sparse factorization. Use SYSMATRIX_MINPIVOTID to retrieve the user equation number that experienced the largest relative drop in its value during factorization.

The matrix determinant is represented by a double precision mantissa times ten raised to an exponent. The exponent is returned using the type SYSMATRIX_EXPONENT. Query for the mantissa using vfs_SysMatrixGetDouble() with type SYSMATRIX_MANTISSA.

Query for the status of the factorization using SYSMATRIX_FACTORIZED. If the function vfs_SysMatrixFactor() fails due to a singular matrix, query for this condition using SYSMATRIX_SINGULAR. If the matrix is singular, the degree of freedom of the singular pivot may be returned using SYSMATRIX_ZEROPIVOT.

Query for the number of iterations with type SYSMATRIX_NUMITER. During an iterative solution, the value of this parameter is updated at every iteration, making this an ideal parameter to track the iteration procedure with a monitor function.

Use a monitor function to retrieve updated parameters during computations.

During vfs_SysMatrixFactor(), if a sparse solver is being used, SYSMATRIX_FACTOR_EQN reflects the monotonically increasing column number being factorized. Likewise, SYSMATRIX_FACTOR_UEQN returns the equation number, in user space, for the current column. These two numbers differ because of reordering that is done for efficiency purposes. If SYSMATRIX_AUTOSPC has been set, then SYSMATRIX_MOD_DIAG returns SYS_ON if the singularity removal algorithm was invoked for this equation.

During vfs_SysMatrixEigen() use SYSMATRIX_EIGENPHASE to determine the phase in process. You can also determine that vfs_SysMatrixEigen() has been called by querying for SYSMATRIX_PHASE and verifying that its value is SYSMATRIX_PHASE_EIGENSOLUTION. If that’s the case, all eigensolver-related queries are available. Depending on the phase, certain parameters may be available. SYSMATRIX_SUBSPACE_SIZE is the size of the subspace used by the eigensolver; SYSMATRIX_EIGENCONV is the number of eigenvalues converged up to the current iteration; SYSMATRIX_EIGEN_STEP is the current subspace iteration number.

The out-of-core solver requires auxiliary memory in the form of a stack. Use SYSMATRIX_STACK_SIZE to query for the number of double precision works used in the stack.

Use SYSMATRIX_PHASE to determine the phase SysMatrix is currently executing. Only two phases are meaningful: SYSMATRIX_PHASE_EIGENSOLUTION and SYSMATRIX_PHASE_SOLVE. That’s because these two phases provide the user with additional information.

During vfs_SysMatrixSolve() - which can be determined by verifying that SYSMATRIX_PHASE takes the value SYSMATRIX_PHASE_SOLVE - you can query for additional iteration parameters depending on the iteration iteration status - SYSMATRIX_ITERSTATUS. If the status is set to SYSMATRIX_ITERSTATUS_INIT then the iteration is about to start; if set to SYSMATRIX_ITERSTATUS_ITER then SysMatrix is iterating to find a solution. If set to SYSMATRIX_ITERSTATUS_TERM then the iteration has terminated. During the iterative solution you also query for SYSMATRIX_NUMITER, the number of iterations. You can also use vfs_SysMatrixGetDouble() to query for specific convergence norms: SYSMATRIX_ITER_FCONV, for residual convergence; SYSMATRIX_ITER_UCONV, for solution convergence; and SYSMATRIX_ITER_ECONV, for energy convergence. If the GMRES solver is used only SYSMATRIX_ITER_FCONV is available.

Errors

SYS_ERROR_ENUM is generated if an improper type is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Type of integer information to query 

    x=SYSMATRIX_EIGENCONV     Number of converged eigenvalues
 =SYSMATRIX_EIGENPHASE    Phase in an eigen solution
 =SYSMATRIX_EIGENSTEP     The subspace iteration number
 =SYSMATRIX_EXPONENT      Exponent base ten of matrix determinant
 =SYSMATRIX_FACTORIZED    Factorized flag
 =SYSMATRIX_FACTOR_EQN    Current factor column from 1 to neq
 =SYSMATRIX_FACTOR_UEQN   Current user-numbered factor equation
 =SYSMATRIX_ITERSTATUS    Status of the iterative solution
 =SYSMATRIX_MAXSUPERSIZE  Size of largest supernode
 =SYSMATRIX_MINPIVOTID    User-numbered minimum pivot equation
 =SYSMATRIX_MOD_DIAG      Toggle for modified factored pivot
 =SYSMATRIX_NUMITER       Number of iterative solver iterations
 =SYSMATRIX_NUMNEG        Number negative roots
 =SYSMATRIX_PHASE         Phase object is in
 =SYSMATRIX_NUM_REQEIGEN  Number of requested eigenvalues
 =SYSMATRIX_SINGULAR      Matrix singular flag
 =SYSMATRIX_STACK_SIZE    Stack size for out-of-core solver
 =SYSMATRIX_SUBSPACE_SIZE Computed eigenvalue subspace size
 =SYSMATRIX_ZEROPIVOT     Zero pivot degree of freedom

  • param – [out] Returned integer information

void vfs_SysMatrixGetLong(vfs_SysMatrix *p, Vint type, Vlong param[])

get long matrix information

Query for long information associated with the system matrix.

The number of non-zero entries in a sparse matrix can be retrieved with SYSMATRIX_ASSEM_NONZ.

The query SYSMATRIX_FACTOR_NONZ requires the system matrix to have been processed using vfs_SysMatrixProcess().

Errors

  • SYS_ERROR_ENUM is generated if an improper type is specified.

  • SYS_ERROR_OPERATION is generated if this function is called with

  • SYSMATRIX_FACTOR_NONZ and the matrix has not been processed using vfs_SysMatrixProcess().

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Type of integer information to query 

    x=SYSMATRIX_ASSEM_NONZ    Number non zeros in assembled matrix
 =SYSMATRIX_FACTOR_NONZ   Number Non zeros in factored matrix

  • param – [out] Returned integer information

void vfs_SysMatrixGetDouble(vfs_SysMatrix *p, Vint type, Vdouble param[])

get double precision matrix information

Query for double precision information associated with the system matrix. If the value has not been computed it is returned as zero.

The double precision queries SYSMATRIX_FACTOR_OPS and SYSMATRIX_SOLVE_OPS require the system matrix to have been processed using vfs_SysMatrixProcess(). The double precision query SYSMATRIX_MANTISSA requires the system matrix to have been factorized using vfs_SysMatrixFactor(). The matrix determinant is represented by a double precision mantissa times ten raised to an exponent. The mantissa is returned using the type SYSMATRIX_MANTISSA. Query for the exponent using vfs_SysMatrixGetInteger() with type SYSMATRIX_EXPONENT.

The double precision queries SYSMATRIX_ITER_FCONV, SYSMATRIX_ITER_UCONV and SYSMATRIX_ITER_ECONV return the normalized residual, solution and energy convergence numbers. These numbers are updated at every iteration, and can be queried in a monitor function. They can be directly compared with the tolerances specified with vfs_SysMatrixSetParamd() using the types SYSMATRIX_ITER_FTOL, SYSMATRIX_ITER_UTOL and SYSMATRIX_ITER_ETOL respectively.

The queries SYSMATRIX_OLD_DIAG and SYSMATRIX_NEW_DIAG return the pivot values before and after factorization. If SYSMATRIX_AUTOSPC is set to SYS_ON, then SYSMATRIX_MOD_DIAG gives the value of the algorithmically modified pivot. All three values are best used in a monitor function along with the integer parameters SYSMATRIX_FACTOR_EQN and SYSMATRIX_FACTOR_UEQN.

Use SYSMATRIX_EIGENSHIFT to retrieve the currently used eigenvalue shift.

Use SYSMATRIX_EIGENVAL and SYSMATRIX_EIGENERR to return an array dimensioned the subspace size with the value of each eigenvalue in the subspace, converged or not, within the required limits or not, and its error, respectively. The subspace size can be obtained by using vfs_SysMatrixGetInteger() and using the parameter SYSMATRIX_SUBSPACE_SIZE.

After factorization, use SYSMATRIX_MINPIVOTVAL to retrieve the minimum relative pivot value drop.

Errors

SYS_ERROR_ENUM is generated if an improper type is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Type of integer information to query 

    x=SYSMATRIX_ITER_FCONV     Iterative solver residual error
 =SYSMATRIX_ITER_UCONV     Iterative solver solution error
 =SYSMATRIX_ITER_ECONV     Iterative solver energy error
 =SYSMATRIX_EIGENSHIFT     The currently used eigenvalue shift
 =SYSMATRIX_EIGENVAL       Eigenvalues in subspace
 =SYSMATRIX_EIGENERR       Error in each eigenvalue in subspace
 =SYSMATRIX_FACTOR_OPS     Number of FLOPS for factor
 =SYSMATRIX_MANTISSA       Mantissa of matrix determinant
 =SYSMATRIX_MINPIVOTVAL    The smallest pivot reduction
 =SYSMATRIX_MOD_DIAG       The pivot's algorithmically updated value
 =SYSMATRIX_NEW_DIAG       The factorized value of the diagonal pivot
 =SYSMATRIX_OLD_DIAG       The original value of the diagonal pivot
 =SYSMATRIX_SOLVE_OPS      Number of FLOPS for solve

  • param – [out] Returned double precision information

void vfs_SysMatrixSetObject(vfs_SysMatrix *p, Vint objecttype, Vobject *object)

set attribute object

Set a pointer to an attribute object. The DofTab object must be registered for symmetrix and non-symmetric sparse and iterative matrices. It must also be entered if any MPC’s are to be applied. Note that the DofTab object can not be set on a diagonal matrix.

Errors

  • SYS_ERROR_OBJECTTYPE is generated if an improper objecttype is specified.

  • SYS_ERROR_OPERATION is generated if the matrix type is diagonal.

Parameters:

  • p – Pointer to SysMatrix object.

  • objecttype – The object type identifier 

    x=VFS_DOFTAB             \docref2{DofTab} object

  • object – Pointer to the object to be set.

void vfs_SysMatrixSetFunction(vfs_SysMatrix *p, Vint functype, VSysMatrixfunc *function, Vobject *object)

set a function pointer

Set a pointer to a function. By default, both the monitor and the matrix product callbacks are NULL. A callback is not invoked if it is NULL.

The method vfs_SysMatrixAbort() can be called from within the monitoring function if necessary.

The available data SysMatrix can be queried for depends on the function called. Functions for which a monitor function is called are vfs_SysMatrixFactor(), vfs_SysMatrixSolve(), and vfs_SysMatrixEigen().

A monitor function can be used with vfs_SysMatrixFactor() if vfs_SysMatrixDef() is defined with SYSMATRIX_SYMM_SPARSE, SYSMATRIX_NSYMM_SPARSE, SYSMATRIX_SYMM_FULL, or SYSMATRIX_SYMM_FULL, SYSMATRIX_TECH is set to SYSMATRIX_TECH_SPARSE (i.e., a direct, sparse factorization is being performed). The function vfs_SysMatrixGetInteger() can be called to query for the parameter SYSMATRIX_FACTOR_EQN, the column number currently being factored. This parameter varies from 1 to the number of equations and can be used to measure progress in the factorization. The user equation can also be retrieved via the parameter SYSMATRIX_FACTOR_UEQN. You can also query for three different double precision parameters with vfs_SysMatrixGetDouble(): SYSMATRIX_OLD_DIAG, SYSMATRIX_NEW_DIAG, and SYSMATRIX_UPD_DIAG. The first parameter gives you the original value of the diagonal for the equation being processed; the second gives its value after factorization. If SYSMATRIX_AUTOSPC is set to SYS_ON, the parameter SYSMATRIX_UPD_DIAG gives the value of the updated pivot. In addition, the integer parameter SYSMATRIX_MOD_DIAG flags whether the pivot has been modified by the presence of SYSMATRIX_AUTOSPC.

During factorization, the frequency with which the monitor function is called can be controlled with the integer parameter SYSMATRIX_MONITOR_FREQ, which defaults to 1000.

When used along with vfs_SysMatrixEigen(), use vfs_SysMatrixGetInteger() with the parameter SYSMATRIX_EIGENPHASE to determine the current phase in the eigensolver. The following phases are available:

Phase

Description

SYSMATRIX_EIGENPHASE_FACTOR

Matrix being factored

SYSMATRIX_EIGENPHASE_SETUP

Factorization complete

SYSMATRIX_EIGENPHASE_ITER

Subspace iteration

During factorization, SYSMATRIX_FACTOR_EQN can be queried for. During the setup phase, vfs_SysMatrixGetInteger() can be used to query for the subspace size with SYSMATRIX_SUBSPACE_SIZE. During the subspace iteration process, vfs_SysMatrixGetInteger() can be used to query for the iteration number with SYSMATRIX_EIGENSTEP, and for the number of converged eigenvalues at the step, with SYSMATRIX_EIGENCONV. At any time during vfs_SysMatrixEigen() you can retrieve the number of eigenvalues requested using vfs_SysMatrixGetInteger() with the parameter SYSMATRIX_NUM_REQEIGEN.

If SYSMATRIX_FUN_PROD is used, the method vfs_SysMatrixGetSysVector() must be called from within the product function; the product is constructed as sysvectorb = A sysvectora. If there are MPCs set, these will be taken care of internally by SysMatrix.

Errors

SYS_ERROR_ENUM is generated if an improper functype is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • functype – The function type identifier 

    x=SYSMATRIX_FUN_MONITOR             Monitoring function
 =SYSMATRIX_FUN_PROD             External function that multiplies the system matrix by a vector

  • function – Pointer to the function to be called

  • object – Pointer to a user object to be passed to the monitoring function.

void vfs_SysMatrixAbort(vfs_SysMatrix *p)

abort a sparse direct factorization

Use vfs_SysMatrixAbort() to abort a sparse, direct factorization in progress. This function is usually called from within a user-defined monitoring function specified with vfs_SysMatrixSetFunction(). The monitoring function is called at regular intervals during the factorization, thereby allowing the user to issue the abort command if necessary.

Parameters:

p – Pointer to SysMatrix object.

void vfs_SysMatrixProcess(vfs_SysMatrix *p)

process before factorization

Perform pre-processing operations on the system matrix required before factorization. If a direct solver option is used, the equations are reordered and a symbolic factorization is performed. If an iterative solver has been specified, this function simply returns.

Errors

SYS_ERROR_OPERATION is generated if the matrix has not been assembled.

Parameters:

p – Pointer to SysMatrix object.

void vfs_SysMatrixReact(vfs_SysMatrix *p, vfs_SysVector *sysvectorx, vfs_SysVector *sysvectory)

compute reaction forces from applied displacements

Compute the reaction forces associated with the degrees of freedom on which a temporary restraint has been set with vfs_SysMatrixSetTemp(). The applied displacements are specified in the system vector sysvectorx and the reactions are stored in the sysvectory.

Parameters:
void vfs_SysMatrixGetSysVector(vfs_SysMatrix *p, Vint type, vfs_SysVector **v)

retrieve internal SysVector objects

Errors

SYS_ERROR_ENUM is generated if an improper type is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • type – Type of SysVector to retrieve 

    x=SYSMATRIX_VEC_PRODA    Vector a
 =SYSMATRIX_VEC_PRODB    Vector b

  • v – [out] Pointer to pointer of the retreived SysVector Retrieve an internal SysVector object. Two vectors are provided: a and b. These vectors should be used if an external function is being used to compute the matrix product b = Aa.

void vfs_SysMatrixGather(vfs_SysMatrix *p, Vint lmcol, Vint *ndofs, Vint lm[], Vdouble s[])

gather non-zero values of a column

Extract the non-zero values of a column of sysmatrix. The column number is in the interval 1 <= lmcol <= number of global dof. The lm and s arrays should be sized to possibly hold the number of global dof defined for the matrix by vfs_SysMatrixDef().

Errors

SYS_ERROR_VALUE is generated if an improper lmcol is specified.

Parameters:

  • p – Pointer to SysMatrix object.

  • lmcol – Column number to extract.

  • ndofs – [out] Number of non-zero degrees of freedom

  • lm – [out] Vector of global dof numbers of non-zero values

  • s – [out] Gathered values

void vfs_SysMatrixPreProcess(vfs_SysMatrix *p)

establish matrix structure

Use this function to establish the matrix structure. It must be called after the DofTab object has been set in SysMatrix. It must be called before any element matrices can be assembled.

Errors

SYS_ERROR_MEMORY is generated if there isn’t enough memory to store the matrix.

Parameters:

p – Pointer to SysMatrix object.

void vfs_SysMatrixSetNumThreads(vfs_SysMatrix *p, Vint num)

set number of threads

Set the number of threads used for parallel execution of SysMatrix functions. If num is 1, then SysMatrix functions are run serially. If num is greater than or equal to 2 then num threads are used for parallel execution.

Errors

SYS_ERROR_VALUE is generated if num is less than 0.

Parameters:

  • p – Pointer to SysMatrix object.

  • num – Number of threads 

    x=0                   Default, serial execution.
 >=1                  Number of threads for parallel execution

void vfs_SysMatrixSetNumThreadsAssem(vfs_SysMatrix *p, Vint nthreads)

set number of threads for assembly

Set the maximum number of threads to be used during parallel assembly.

Errors

SYS_ERROR_VALUE is generated if nthreads is less than one.

Parameters:

  • p – Pointer to SysMatrix object.

  • nthreads – Maximum number of threads to be used.

void vfs_SysMatrixSetInitSol(vfs_SysMatrix *p, vfs_SysVector *initsol)

set initial solution vector

Set an initial guess solution vector to be used with the iterative solver.

Errors

SYS_ERROR_OPERATION is generated if the number of global dof in system matrix and system vector are not equal.

Parameters:

  • p – Pointer to SysMatrix object.

  • initsol – Pointer to SysVector object initial guess solution

void vfs_SysMatrixPre(vfs_SysMatrix *p, Vint complextype)

specify real or complex matrix

Use vfs_SysMatrixPre() to specify whether the matrix is real or complex. For real matrices complextype takes the value SYS_DOUBLE. For complex matrices, complextype takes the value SYS_DOUBLECOMPLEX.

This call must be made before vfs_SysMatrixDef() is called. Defaults to real.

Errors

SYS_ERROR_ENUM is issued if complextype is invalid.

Parameters:

  • p – Pointer to SysMatrix object.

  • complextype – Enumerated type for real or complex matrix 

    x=SYS_DOUBLE              Real matrix
 =SYSMATRIX_DOUBLECOMPLEX Complex matrix

void vfs_SysMatrixSetComplexMode(vfs_SysMatrix *p, Vint complexmode)

set complex mode

Specify floating point data mode.

The methods vfs_SysMatrixAssem(), vfs_SysMatrixAssemDiag(), vfs_SysMatrixAssemFull(), and all take arrays of floating point numbers as lower triangular, full, or diagonal element matrices. In addition, vfs_SysMatrixColumn() returns an array of floating point numbers for the column values. The values entered or returned can be declared as SYS_COMPLEX_REAL, SYS_COMPLEX_IMAGINARY, or SYS_COMPLEX_REALIMAGINARY. The data in the argument then takes the specified mode.

This is useful if a real element stiffness matrix is being assembled. The mode SYS_COMPLEX_REAL indicates that only one value is being entered for each a_ij in the array, and this value is assembled into the matrix. In this case, the imaginary entries in the matrix remain unchanged. Likewise, if an imaginary damping matrix is being assembled, the mode SYS_COMPLEX_IMAGINARY can be specified so that only the imaginary parts of SysMatrix are updated. If mode is set to SYS_COMPLEX_REALIMAGINARY then two consecutive numbers in the array stand for the real and imaginary parts of the data. Defaults to SYS_COMPLEX_REAL.

Parameters:

  • p – Pointer to SysMatrix object.

  • complexmode – The type of data being passed through the Vdouble[] argument list. 

    x=SYS_COMPLEX_REAL          Data is real
 =SYS_COMPLEX_IMAGINARY     Data is imaginary
 =SYS_COMPLEX_REALIMAGINARY Data is complex

void vfs_SysMatrixSetThreadComplexMode(vfs_SysMatrix *p, Vint ithread, Vint complexmode)

set complex mode in thread

Specify floating point data mode.

The methods vfs_SysMatrixThreadAssem(), vfs_SysMatrixThreadAssemDiag(), and vfs_SysMatrixThreadAssemFull() all take arrays of floating point numbers as lower triangular, full, or diagonal element matrices. In addition, vfs_SysMatrixColumn() returns an array of floating point numbers for the column values. The values entered or returned can be declared as SYS_COMPLEX_REAL, SYS_COMPLEX_IMAGINARY, or SYS_COMPLEX_REALIMAGINARY. The data in the argument then takes the specified mode.

This is useful if a real element stiffness matrix is being assembled. The mode SYS_COMPLEX_REAL indicates that only one value is being entered for each a_ij in the array, and this value is assembled into the matrix. In this case, the imaginary entries in the matrix remain unchanged. Likewise, if an imaginary damping matrix is being assembled, the mode SYS_COMPLEX_IMAGINARY can be specified so that only the imaginary parts of SysMatrix are updated. If mode is set to SYS_COMPLEX_REALIMAGINARY then two consecutive numbers in the array stand for the real and imaginary parts of the data. Defaults to SYS_COMPLEX_REAL.

Errors

  • SYS_ERROR_ENUM is generated if an improper complexmode is specified.

  • SYS_ERROR_OPERATION is generated if the matrix is not defined as complex.

  • SYS_ERROR_OPERATION is generated if vfs_SysMatrixSetNumThreadsAssem() has not been previously called.

  • SYS_ERROR_VALUE is generated if ithread is less than 1 or greater than the maximum number of threads specified in vfs_SysMatrixSetNumThreadsAssem().

Parameters:

  • p – Pointer to SysMatrix object.

  • ithread – Thread number starting from 1.

  • complexmode – The type of data being passed through the Vdouble[] argument list. 

    x=SYS_COMPLEX_REAL          Data is real
 =SYS_COMPLEX_IMAGINARY     Data is imaginary
 =SYS_COMPLEX_REALIMAGINARY Data is complex