.x.
|
|
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
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 Mu = w. 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
