# 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 object. Memory is allocated for the object private data and the pointer to the object is returned.

Destroy an instance of a object using

```void vfs_SysVectorEnd (vfs_SysVector *sysvector)
```

Return the current value of a object error flag using

```Vint vfs_SysVectorError (vfs_SysVector *sysvector)
```

Make a copy of a 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 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 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 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 object and moved to another object. If each random vector were to be generated in separate objects they would all be random but identical.

Parameters

p – Pointer to 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 object.

• nedofs – Number of element degrees of freedom

• lm – Vector of element global dof numbers

• v – Element vector

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 object.

• 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 object.

• ndofs – Number of degrees of freedom

• lm – Vector of global dof numbers

• v[out] Gathered values

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 object.

• 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 = v
a_lm = v
...
a_lm[ndofs-1] = v[ndofs-1]
```
Errors

`SYS_ERROR_VALUE` is generated if an improper ndofs is specified.

Parameters
• p – Pointer to object.

• ndofs – Number of degrees of freedom

• lm – Vector of global dof numbers

• v – Scattered values

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 object.

• 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 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 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 object.

• sysvectorx – Pointer to 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 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 object.

• sysvectorx – Pointer to 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 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 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, sysvectorx and sysvectory to sysvector . The operation replaces any previous values in with the sum of sysvectorx and s times sysvectory. A single 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 object.

• sysvectorx – Pointer to object.

• s – Scale factor

• sysvectory – Pointer to 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
a_2 = sysvector * sysvectort
...
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 object.

• nvec – Number of system vectors 0 < nvec

• sysvectort – Poinsysvectorter to array of nvec objects

• sysvectorq – Pointer to reduced 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 + a_2 * sysvectort + ... + 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 object.

• nvec – Number of system vectors 0 < nvec

• sysvectort – Pointer to array of nvec objects

• sysvectorq – Pointer to expanded 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
• p – Pointer to object.

• sysvectorx – Pointer to object x.

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
• p – Pointer to object.

• sysvectorx – Pointer to object.

• sysvectory – Pointer to object.

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 object.

• s – Scale factor

• sysvectorx – Pointer to object.

• sysvectory – Pointer to 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
• p – Pointer to object.

• sysvectorx – Pointer to object x.

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 object.

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

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
...
```

Parameters
• p – Pointer to object.

• sysvectorx – Pointer to 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.

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 = x_1
...
```

Parameters
• p – Pointer to object.

• sysvectorx – Pointer to 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.

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
a_2 = x_2 + s * y_mask
...
a_[n-1] = x_[n-1] + s * y_mask[n-1]
```

Parameters
• p – Pointer to object.

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

• s – Scaling parameter.

• sysvectory – Pointer to 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 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 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
```

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 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
```

Set the number of threads used for parallel execution of functions. If num is 1, then 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 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;

/* 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);

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 object. Memory is allocated for the object private data and the pointer to the object is returned.

Destroy an instance of a object using

```void vfs_SysMatrixEnd (vfs_SysMatrix *sysmatrix)
```

Return the current value of a object error flag using

```Vint vfs_SysMatrixError (vfs_SysMatrix *sysmatrix)
```

Returns

The function returns a pointer to the newly created 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 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 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 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 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 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 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()`. 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 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 object.

• nedofs – Number of element degrees of freedom

• lm – Vector of element global dof numbers

• s – Element matrix

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 object.

• 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 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 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 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
• p – Pointer to object.

• sysvectorx – Pointer to object.

• sysvectory – Pointer to object.

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 object.

• nvec – Number of system vectors 0 < nvec

• sysvectorx – Pointer to array of nvec objects

• sysvectory – Pointer to array of nvec 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 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 object.

• lmcol – Column number to extract.

• sysvector – Pointer to 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
• p – Pointer to object.

• sysvectorx – Pointer to object.

• sysvectory – Pointer to object.

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
• p – Pointer to object.

• sysvector – Pointer to object

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 = [ a_00 a_01 ... a_0N ]
sysvector = [ 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 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 objects

• sysmatrixq – Pointer to 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 object.

• nvec – number of vectors X and Y.

• sysvectorx – Pointer to vector of objects X

• sysvectory – Pointer to vector of 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 object.

• sysmatrixb – Pointer to object for matrix B.

• sysvectore – Pointer to object of eigenvalues

• sysvectorx – Pointer to array of 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 object.

• sysmatrixb – Pointer to 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 object.

• ieigen – Eigen mode to retrieve

• sysvector – Pointer to 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 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 object.

• s – number of sums

• sysmatrixb – Pointer to 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 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 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 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 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 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 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 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 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 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 .

Errors

`SYS_ERROR_ENUM` is generated if an improper functype is specified.

Parameters
• p – Pointer to 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 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 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
• p – Pointer to object.

• sysvectorx – Pointer to object.

• sysvectory[out] Pointer to object.

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 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 Retrieve an internal 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 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 object has been set in . 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 object.

Set the number of threads used for parallel execution of functions. If num is 1, then 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 object.

• num – Number of threads

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

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 object.

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 object.

• initsol – Pointer to 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 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 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 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
```

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 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 object.

```x=SYS_COMPLEX_REAL          Data is real