Goal

Obtain a solution to a boundary value problem for the thermal equation, with thermal coefficients that depend on the solution.

Solution

Use a fixed-point iteration approach [Amos10], utilizing Intel MKL PARDISO for solving linear problems on each external iteration.

  1. Set up the matrix structure in CSR format.

  2. Perform fixed-point iteration until the residual norm becomes lower than the tolerance.

    1. Use the pardiso routine to solve the linearized system for the current iteration.

    2. Set the solution of the system to the next approximation of the main equation using the dcopy routine.

    3. Based on the new approximation, calculate the new elements of the matrix.

    4. Calculate the residual of the current solution using the mkl_cspblas_dcsrgemv routine.

    5. Calculate the norm of the residual using the dnrm2 routine and compare it with the tolerance.

  3. Free the internal memory of the solver.

Source code: see the sparse folder in the samples archive available at http://software.intel.com/en-us/mkl_cookbook_samples.

Finding an approximate solution using Intel MKL PARDISO, Sparse BLAS, and BLAS

      CONSTRUCT_MATRIX_STRUCTURE (nx, ny, nz, &ia, &ja, &a, &error);
      CONSTRUCT_MATRIX_VALUES (nx, ny, nz, us, ia, ja, a, &error);
      DO WHILE res > tolerance
        phase = 13;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, 
          &nrhs, iparm, &msglvl, f, u, &error );      
        DCOPY (&n, u, &one, u_next, &one);
        CONSTRUCT_MATRIX_VALUES (nx, ny, nz, u_next, ia, ja, a, &error);
        MKL_CSPBLAS_DCSRGEMV (uplo, &n, a, ia, ja, u, temp );
        DAXPY (&n, &minus_one, f, &one, temp, &one);
        res = DNRM2 (&n, temp, &one);
      END DO
      phase = -1;
      PARDISO ( pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum,
        &nrhs, iparm, &msglvl, f, u, &error );

Routines Used

Task

Routine

Description

Solve the linearized system on the current iteration; free internal memory of the solver.

PARDISO

Calculates the solution of a set of sparse linear equations with multiple right-hand sides.

Set the solution found as the next approximation of the main equation.

DCOPY

Copies vector to another vector.

Calculate the residual of the current nonlinear iteration.

MKL_CSPBLAS_DCSRGEMV

Computes matrix - vector product of a sparse general matrix stored in the CSR format (3-array variation) with zero-based indexing.

DAXPY

Computes a vector-scalar product and adds the result to a vector.

Calculate the norm of the residual to compare it with stopping criteria.

DNRM2

Computes the Euclidean norm of a vector.

Discussion

The stationary nonlinear heat equation can be described as a boundary value problem for a nonlinear partial differential equation:



Where the domain D is assumed to be a cube: , and is an unknown function of temperature.

For the purpose of demonstration, the problem is restricted to linear dependence of the thermal coefficient on the solution:



To obtain a numerical solution, an equidistant grid with grid step h in the domain D is chosen, and the partial differential equation is approximated using finite differences. This procedure [Smith86] yields a system of nonlinear algebraic equations:









Each equation ties together the value of the unknown grid function u and the value of the respective right hand side at seven grid points. The left hand sides of the equations can be represented as linear combinations of the grid function values with coefficients which depend on the solution itself. Introducing a matrix composed of these coefficients, the equations can be rewritten in vector-matrix form:



Since the coefficient matrix A is sparse (it has only seven nonzero elements in each row), it is suitable to store it in a CSR-format array (see Sparse Matrix Storage Formats in the Intel Math Kernel Library Reference Manual ), and use the PARDISO* solver for solving it using this iterative algorithm:

  1. Set u to initial value u 0.

  2. Calculate residual r = A(u)u - g.

  3. Do while ||r|| < tolerance:

    1. Solve system A(u)w = g for w.

    2. Set u = w.

    3. Calculate residual r = A(u)u - g.

Для получения подробной информации о возможностях оптимизации компилятора обратитесь к нашему Уведомлению об оптимизации.