Computes matrix-matrix or matrix-vector product.

Operator Return Value

Matrix containing the product of A and B.   (Output)

Required Operands

A — Left operand matrix or vector. This is an array of rank 1, 2, or 3. It may be real, double, complex, double complex, or one of the computational sparse matrix derived types, ?_hbc_sparse. (Input)
Note that A and B cannot both be ?_hbc_sparse.

B — Right operand matrix or vector. This is an array of rank 1, 2, or 3. It may be real, double, complex, double complex, or one of the computational sparse matrix derived types, ?_hbc_sparse. (Input)
Note that A and B cannot both be ?_hbc_sparse.

If A has rank one, B must have rank two.
If B has rank one, A must have rank two.
If A has rank three, B must have rank three.
If B has rank three, A must have rank three.

FORTRAN 90 Interface

A .x. B

Description

Computes the product of  matrix or vector A and matrix or vector B. The results are in a precision and data type that ascends to the most accurate or complex operand.

Rank three operation is defined as follows:

do i = 1, min(size(A,3), size(B,3))
   X(:,:,i) =  A(:,:,i) .x. B(:,:,i)
end do

 

.x. can be used with either dense or sparse matrices. It is MPI capable for dense matrices only.

Examples

Dense Matrix Example  (operator_ex03.f90)

 

  use linear_operators

      implicit none

 

! This is the equivalent of Example 3 for LIN_SOL_GEN using operators.

      integer, parameter :: n=32

      real(kind(1e0)) :: one=1e0, zero=0e0, A(n,n), b(n), x(n)

      real(kind(1e0)) change_new, change_old

      real(kind(1d0)) :: d_zero=0d0, c(n), d(n,n), y(n)

 

! Generate a random matrix and right-hand side.

      A = rand(A); b= rand(b)

 

! Save double precision copies of the matrix and right-hand side.

      D = A

      c = b

! Compute single precision inverse to compute the iterative refinement.

      A = .i. A

 

! Start solution at zero.  Update it to an accurate solution

! with each iteration.

      y = d_zero

      change_old = huge(one)

 

      iterative_refinement: do

! Compute the residual with higher accuracy than the data.

         b = c - (D .x. y)

 

! Compute the update in single precision.

         x = A .x. b

         y = x + y

         change_new = norm(x)

 

! Exit when changes are no longer decreasing.

         if (change_new >= change_old) exit iterative_refinement

         change_old = change_new

      end do iterative_refinement

 

      write (*,*) 'Example 3 for LIN_SOL_GEN (operators) is correct.'

      end

Sparse Matrix Example

Consider the one-dimensional Dirichlet problem

Using a standard approach to solving this involves approximating the second derivative operator with central divided differences

This leads to the sparse linear algebraic system.  The definitions for these terms are implied in the following Fortran program.

 

Subroutine document_ex1

! Illustrate a 1D Poisson equation with Dirichlet boundary conditions.

! This module defines the structures and overloaded assignment code.

      Use linear_operators

      Implicit None

!

      Integer :: I

      Integer, Parameter :: N = 1000

      Real (Kind(1.d0)) :: f, h, r, w (N), a = 0.d0, b = 1.d0, &
      u_a = 0.d0, u_b = 1.d0, u (N)

      Type (d_sparse) M

      Type (d_hbc_sparse) K

      External f

! Define the difference used.

      h = (b-a) / (N-1)

      r = 1.d0 / h ** 2

! Fill in the matrix entries.

! Isolated equation for the left boundary condition.

      M = d_entry (1, 1, r)

      Do I = 2, N - 1

         M = d_entry (I, I-1, r)

         M = d_entry (I, I,-2*r)

         M = d_entry (I, I+1, r)

      End Do

      

! Isolated equation for the right boundary condition. 

      M = d_entry (N, N, r)

      

! Fill in the right-hand side (a dense vector).

      Do I = 2, N - 1

         w (I) = f (a+(I-1)*h)

      End Do

! Insert the known end conditions.  These should be satisfied

! almost exactly, up to rounding errors.

      w (1) = u_a * r

      w (N) = u_b * r

! Ready to solve …

! Conversion to Harwell-Boeing format using overloaded assignment

 

      K = M

! Solve the system using an IMSL defined operator.

 

      u = K .ix. w  

! The parentheses are needed because of precedence rules.

! Compute residuals and overwrite w(:) with these values.

             

      w = w - (K .x. u)

End Subroutine

!

Function f (x)

      Real (Kind(1.d0)) :: f, x

! Define a hat function, peaked at x=0.5.

      If (x <= 0.5d0) Then

         f = x

      Else

         f = 1.d0 - x

      End If

End Function

Parallel Example  (parallel_ex03.f90)

This example shows the box data type used while obtaining an accurate solution of several systems.  Important in this example is the fact that only the root will achieve convergence, which controls program flow out of the loop.  Therefore the nodes must share the root's view of convergence, and that is the reason for the broadcast of the update from root to the nodes.  Note that when writing an explicit call to an MPI routine there must be the line INCLUDE ‘mpif.h', placed just after the IMPLICIT NONE statement.  Any number of nodes can be used.

 

      use linear_operators

      use mpi_setup_int

 

      implicit none

      INCLUDE 'mpif.h'

 

! This is the equivalent of Parallel Example 3 for .i. and iterative

! refinement with box date types, operators and functions.

      integer, parameter :: n=32, nr=4

      integer IERROR

      real(kind(1e0)) :: one=1e0, zero=0e0

      real(kind(1e0)) :: A(n,n,nr), b(n,1,nr), x(n,1,nr)

      real(kind(1e0)) change_old(nr), change_new(nr)

      real(kind(1d0)) :: d_zero=0d0, c(n,1,nr), D(n,n,nr), y(n,1,nr)

 

! Setup for MPI.

      MP_NPROCS=MP_SETUP()

    

! Generate a random matrix and right-hand side.

      A = rand(A); b= rand(b)

 

! Save double precision copies of the matrix and right-hand side.

      D = A

      c = b

 

! Get single precision inverse to compute the iterative refinement.

      A = .i. A

 

! Start solution at zero.  Update it to a more accurate solution

! with each iteration.

      y = d_zero

      change_old = huge(one)

 

      ITERATIVE_REFINEMENT: DO

 

! Compute the residual with higher accuracy than the data.

         b = c - (D .x. y)

 

! Compute the update in single precision.

         x = A .x. b

         y = x + y

         change_new = norm(x)

 

! All processors must share the root's test of convergence.

         CALL MPI_BCAST(change_new, nr, MPI_REAL, 0, &

           MP_LIBRARY_WORLD, IERROR)

 

! Exit when changes are no longer decreasing.

         if (ALL(change_new >= change_old)) exit iterative_refinement

         change_old = change_new

      end DO ITERATIVE_REFINEMENT

 

        IF(MP_RANK == 0) write (*,*) 'Parallel Example 3 is correct.'

 

! See to any error messages and quit MPI.

      MP_NPROCS=MP_SETUP('Final')

      end

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260