From GHER

This module used in the Ocean Assimilation Kit implements basic matrix operations which call the corresponding BLAS routines. The operators can be applied to sparse matrices and full matrices.

Operators

A.x.B

A.tx.B

A.xt.B

V.dx.B

A.xd.V

and can either be full matrices (standard Fortran 2-d arrays in single or double precision) or sparse matrices (defined with the Fortran type type(SparseMatrix)).

Functions

The module implements also the functions:

eye(n)

creates an indentify matrix of size n (full matrix)

ones(m,n)

creates m x n matrix with all elements equal to one

zeros(m,n)

creates m x n matrix with all elements equal to zero

inv(A)

computes the inverse of A

det(A)

computes the determinant of A

rand(n,m), randn(n,m)

random matrix (unform distribution or gaussian distribution)

full(S)

converts a sparse matrix to a full matrix

sparse(X)

converts a full matrix to sparse

svd

singular value decomposition

The functions uses an API with is simular to Matlab/Octave whenever possible.

Benchmarks

Machine: Intel Xeon CPU L5420 at 2.50GHz (nic3).

A benchmarks has been conducted for different optimization level calling directly the BLAS routines or using the operators (option -DWITH_OPERATORS). Every test has been repeated 5 times (on a otherwise idle node) and the median cpu time is used.

To run this benchmark you need matoper.F90 and matoper_inc.F90 from the Ocean Assimilation Kit.

Compiler	O0	O0 -DWITH_OPERATORS	O2	O2 -DWITH_OPERATORS	O3	O3 -DWITH_OPERATORS
ifort	2.584	2.593	2.581	2.568	2.582	2.589
pgf90	2.390	2.384	2.372	2.362	2.374	2.363
gfortran	2.585	2.590	2.554	2.552	2.553	2.552

Test code:

program test_matoper
 use matoper
 implicit none
 integer, parameter :: r = 100, m=100000
 real(8) :: lambda(r), UT(r,r), yo(m), Hxf(m), invsqrtR(m), HSf(m,r)

 real(8) :: dummy(1,1)
 integer :: info,i,j

 real(8) :: temp(m,r), ampl(r)
 real(8) :: start, finish

#ifndef WITH_OPERATORS
 real(8) :: tmp(m), tmp2(r), tmp3(r)
#endif 

 do i=1,r
   yo(i) = mod(real(i,8),10._8)
   Hxf(i) = mod(real(i,8),5._8)
   invsqrtR(i) = 6._8
 end do

 do j=1,r
   do i=1,m
     HSf(i,j) =  mod(real(i+j,8),10._8)/10._8 - 0.5_8
     temp(i,j) = invsqrtR(i)*HSf(i,j)     
   end do
 end do


  call gesvd('n','a',temp,lambda,dummy,UT,info)
  lambda = (1+lambda**2)**(-1)

  call cpu_time(start)

  do i=1,100

#ifdef WITH_OPERATORS
    ampl = UT.tx.(lambda.dx.(UT.x.(HSf.tx.(invsqrtR**2*(yo-Hxf)))))
#else
    tmp = invsqrtR**2*(yo-Hxf)
    call DGEMV('T',m,r,1._8,HSf,m,tmp,1,0._8,tmp2,1)
    call DGEMV('N',r,r,1._8,UT,r,tmp2,1,0._8,tmp3,1)
    tmp3 = lambda * tmp3
    call DGEMV('T',r,r,1._8,UT,r,tmp3,1,0._8,ampl,1)
#endif

    ! make loop interation depedent on each other
    yo(1) = yo(1)+ampl(1)/1000._8

  end do

  call cpu_time(finish)

!  write(6,*) 'check sum ',sum(ampl)
  write(6,*) finish-start

end program test_matoper