librsb  1.3.0.1
Functions/Subroutines
fortran.F90 File Reference

Sparse BLAS-based usage: uscr_begin(), usgp(), ussp(), usmv(), ... More...

Functions/Subroutines

subroutine blas_sparse_mod_example (res)
 
subroutine blas_sparse_io_example (res)
 
program main
 

Detailed Description

Sparse BLAS-based usage: uscr_begin(), usgp(), ussp(), usmv(), ...

!
! Copyright (C) 2008-2022 Michele Martone
!
! This file is part of librsb.
!
! librsb is free software; you can redistribute it and/or modify it
! under the terms of the GNU Lesser General Public License as published
! by the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! librsb is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with librsb; see the file COPYING.
! If not, see <http://www.gnu.org/licenses/>.
!
!> @file
!! @brief Sparse BLAS-based usage:
!! uscr_begin(), usgp(), ussp(), usmv(), ...
!! \include fortran.F90
SUBROUTINE blas_sparse_mod_example(res)
USE rsb ! For the second part of the example and RSB_IDX_KIND
USE iso_c_binding
IMPLICIT NONE
INTEGER :: i, j
INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0, res
TYPE(C_PTR),TARGET :: mtxAp = c_null_ptr ! matrix pointer
INTEGER :: A
INTEGER,PARAMETER :: transn = blas_no_trans
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
REAL(KIND=8),parameter :: alpha = 3
! Symmetric (declared via lower triangle) matrix based example, e.g.:
! 1 0
! 1 1
! declaration of VA,IA,JA
!INTEGER,PARAMETER :: nr = 100
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nr = 20
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nc = nr
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nnz = (nr*(nr+1))/2 ! half the square
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 1
INTEGER(KIND=RSB_IDX_KIND) :: nt = 0
INTEGER :: ic, ir
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: IA(nnz) = (/ (((ir), ic=1,ir), ir=1,nr ) /) ! (/1, 2, 2/)
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: JA(nnz) = (/ (((ic), ic=1,ir), ir=1,nr ) /) ! (/1, 1, 2/)
REAL(KIND=8),parameter :: va(nnz) = (/ ((1, ic=1,ir), ir=1,nr ) /) ! (/1, 1, 1/)
REAL(KIND=8) :: x(nc,nrhs) = reshape((/((1), ic=1,nc*nrhs)/),[nc,nrhs]) ! reference x ! (/1, 1/)
REAL(KIND=8),parameter :: cy(nr,nrhs) = reshape((/((alpha+alpha*nr), ir=1,nr*nrhs)/),[nr,nrhs]) ! reference cy after ! (/9, 9/)
REAL(KIND=8) :: y(nr,nrhs) = reshape((/((alpha), ir=1,nr*nrhs)/),[nr,nrhs]) ! y will be overwritten ! (/3, 3/)
! First example part: pure blas_sparse code.
res = 0
CALL duscr_begin(nr,nc,a,res)
IF (res.NE.0) GOTO 9999
CALL ussp(a,blas_lower_symmetric,istat)
IF (istat.NE.0) GOTO 9997
CALL ussp(a,blas_rsb_spmv_autotuning_on,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) print *,"autotuning returned nonzero:", istat &
&," ...did you enable autotuning ?"
!
! First style example
CALL uscr_insert_entries(a,nnz,va,ia,ja,istat)
IF (istat.NE.0) GOTO 9997
CALL uscr_end(a,istat)
IF (istat.NE.0) GOTO 9997
! CALL ussp(A,blas_rsb_duplicates_sum,istat)
! CALL uscr_insert_entries(A,nnz,VA,IA,JA,istat) ! uncomment this to activate add of coefficients to pattern
CALL usgp(a,blas_rsb_spmv_autotuning_on,nt) ! (experimental)
IF (nt.NE.0) print*,"autotuner chose ",nt," threads"
CALL ussp(a,blas_rsb_spmv_autotuning_off,istat) ! (experimental) turns auto-tuning + thread setting off
IF (istat.NE.0) GOTO 9997
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
IF (istat.NE.0) GOTO 9997
!
DO j = 1, nrhs
DO i = 1, nr
IF (y(i,j).NE.cy(i,j)) print *, "first check results are not ok"
IF (y(i,j).NE.cy(i,j)) GOTO 9997
END DO
END DO
!
y(:,:) = alpha ! reset
!
! Second style example
CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) GOTO 9997
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat) ! Equivalent to the above (as long as incx=incy=1).
CALL usmm(blas_colmajor,transn,nrhs,-alpha,a,x,nr,y,nc,istat) ! Subtract the last usmm call contribution.
IF (istat.NE.0) GOTO 9997
!
DO j = 1, nrhs
DO i = 1, nr
IF (y(i,j).NE.cy(i,j)) print *,"second check results are not ok"
IF (y(i,j).NE.cy(i,j)) GOTO 9997
END DO
END DO
!
print *, "check results are ok"
! Second part of the example: access to the rsb.h interface via
! the ISO C Binding interface.
mtxap = rsb_blas_get_mtx(a) ! get pointer to rsb structure (as in the rsb.h API)
IF(nr.LT.5) istat = rsb_file_mtx_save(mtxap,c_null_char) ! write to stdout (only if matrix small enough)
GOTO 9998
9997 res = -1
9998 CONTINUE
CALL usds(a,istat)
IF (istat.NE.0) res = -1
9999 CONTINUE
end SUBROUTINE blas_sparse_mod_example
! Example loading a matrix from file and measuring SPMM.
SUBROUTINE blas_sparse_io_example(res)
USE rsb ! For rsb_blas_file_mtx_load
USE iso_c_binding
IMPLICIT NONE
INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
INTEGER :: j
INTEGER(KIND=RSB_BLAS_IDX_KIND) :: istat = 0
INTEGER :: A
INTEGER,PARAMETER :: transn = blas_no_trans
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incx = 1
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: incy = 1
COMPLEX(KIND=8),PARAMETER :: alpha = 3
INTEGER(KIND=RSB_IDX_KIND) :: nr
INTEGER(KIND=RSB_IDX_KIND) :: nc
INTEGER(KIND=RSB_IDX_KIND) :: nz
INTEGER(KIND=RSB_IDX_KIND) :: st
INTEGER(KIND=RSB_IDX_KIND),PARAMETER :: nrhs = 4
COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: x
COMPLEX(KIND=8), ALLOCATABLE, TARGET, DIMENSION(:,:) :: y
CHARACTER(KIND=C_CHAR,LEN=7),TARGET :: filename = 'pd.mtx'//c_null_char
INTEGER(C_SIGNED_CHAR) :: typecode = rsb_numerical_type_double_complex
REAL(KIND=c_double) :: mvt,mmt,omt
INTEGER(KIND=C_INT),TARGET::IZERO=0
res = 0
a = rsb_blas_file_mtx_load(filename, typecode);
IF (a.EQ.blas_invalid_handle) GOTO 9997
CALL usgp(a,blas_num_rows,nr)
CALL usgp(a,blas_num_cols,nc)
print*,"Read matrix ",filename(1:6)," ",nr,"x",nc,":",nz
CALL usgp(a,blas_general,st)
IF (st .EQ. 1) print*,"Matrix has no symmetry"
IF (st .EQ. 1) print*,"Matrix is upper symmetric"
IF (st .EQ. 1) print*,"Matrix is upper hermitian"
! ...
IF (istat.NE.0) GOTO 9997
WRITE(*,'(a,i0)') "Using NRHS=",nrhs
ALLOCATE( x(nc,nrhs))
ALLOCATE( y(nr,nrhs))
x = 1.0
y = 0.0
mvt = -rsb_time()
DO j = 1, nrhs
CALL usmv(transn,alpha,a,x(:,j),incx,y(:,j),incy,istat)
END DO
IF (istat.NE.0) GOTO 9997
mvt = mvt + rsb_time()
WRITE(*,'(a,e12.4,a)') "Repeated USMV took ",mvt," s"
y = 0.0
mmt = -rsb_time()
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
mmt = mmt + rsb_time()
WRITE(*,'(a,e12.4,a)') "A single USMM took ",mmt," s"
WRITE(*,'(a,g11.4,a)')"USMM-to-USMV speed ratio is is ", mvt/mmt, "x"
print*,"Call auto-tuning routine.."
! Change IZERO value to 1 to get verbose tuning again.
IF (res.NE.0) GOTO 9997
CALL ussp(a,blas_rsb_autotune_next_operation,istat) ! (experimental) turns auto-tuning + thread setting on
IF (istat.NE.0) GOTO 9997
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
print*,"Repeat measurement."
y = 0.0
omt = -rsb_time()
CALL usmm(blas_colmajor,transn,nrhs, alpha,a,x,nr,y,nc,istat)
IF (istat.NE.0) GOTO 9997
omt = omt + rsb_time()
WRITE(*,'(a,e12.4,a)') "Tuned USMM took ",omt," s"
WRITE(*,'(a,g11.4,a)')"Tuned-to-untuned speed ratio is is ",mmt/omt,"x"
GOTO 9998
9997 res = -1
9998 CONTINUE
CALL usds(a,istat)
IF (istat.NE.0) res = -1
end SUBROUTINE blas_sparse_io_example
PROGRAM main
USE rsb, ONLY: rsb_lib_init, rsb_lib_exit, c_ptr, c_null_ptr,&
! USE blas_sparse, only: RSB_BLAS_IDX_KIND ! only if using long indices
USE iso_c_binding
IMPLICIT NONE
INTEGER :: passed = 0, failed = 0
INTEGER(KIND=RSB_IDX_KIND) :: res ! note: same as RSB_BLAS_IDX_KIND
!TYPE(C_PTR),PARAMETER :: EO = RSB_NULL_EXIT_OPTIONS
!TYPE(C_PTR),PARAMETER :: IO = RSB_NULL_INIT_OPTIONS
! Note: using C_NULL_PTR instead of the previous lines because of http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59411
TYPE(C_PTR),PARAMETER :: EO = c_null_ptr
TYPE(C_PTR),PARAMETER :: IO = c_null_ptr
INTEGER(KIND=C_INT),TARGET::IONE=1
res = rsb_lib_init(io)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
res = rsb_lib_exit(eo)
print *, "FAILED:", failed
print *, "PASSED:", passed
IF (failed .GT. 0) THEN
stop 1
END IF
END PROGRAM
auto main() -> int
Definition: assemble.cpp:38
@ blas_no_trans
Definition: blas_sparse.h:32
@ blas_num_nonzeros
Definition: blas_sparse.h:129
@ blas_num_cols
Definition: blas_sparse.h:128
@ blas_num_rows
Definition: blas_sparse.h:127
@ blas_invalid_handle
Definition: blas_sparse.h:133
@ blas_upper_symmetric
Definition: blas_sparse.h:114
@ blas_upper_hermitian
Definition: blas_sparse.h:116
@ blas_general
Definition: blas_sparse.h:107
@ blas_lower_symmetric
Definition: blas_sparse.h:113
@ blas_colmajor
Definition: blas_sparse.h:28
@ blas_rsb_autotune_next_operation
Definition: blas_sparse.h:153
@ blas_rsb_spmv_autotuning_off
Definition: blas_sparse.h:148
@ blas_rsb_spmv_autotuning_on
Definition: blas_sparse.h:147
subroutine blas_sparse_io_example(res)
Definition: fortran.F90:126
subroutine blas_sparse_mod_example(res)
Definition: fortran.F90:27
Definition: rsb_blas_sparse.F90:243
inserts multiple entries
Definition: rsb_blas_sparse.F90:49
multiplication : c <- beta c + alpha A b
Definition: rsb_blas_sparse.F90:126
multiplication : c <- beta c + alpha A b
Definition: rsb_blas_sparse.F90:104
ISO C BINDING interface to rsb_blas_file_mtx_load.
Definition: rsb.F90:864
ISO C BINDING interface to rsb_file_mtx_save.
Definition: rsb.F90:745
ISO C BINDING interface to rsb_lib_exit.
Definition: rsb.F90:107
ISO C BINDING interface to rsb_lib_init.
Definition: rsb.F90:49
ISO C BINDING interface to rsb_lib_set_opt.
Definition: rsb.F90:83
ISO C BINDING interface to rsb_time.
Definition: rsb.F90:851
Definition: rsb_blas_sparse.F90:29
subroutine usgp(A, pname, istat)
Get a matrix property.
Definition: rsb_blas_sparse.F90:301
subroutine ussp(A, pname, istat)
Set a matrix property. Should be called just after creation, before nonzeroes insertion.
Definition: rsb_blas_sparse.F90:324
subroutine uscr_end(A, istat)
Makes an assembled matrix out of a matrix in build state. After this, it is not possible anymore to i...
Definition: rsb_blas_sparse.F90:279
subroutine usds(A, istat)
Destroys a matrix.
Definition: rsb_blas_sparse.F90:259
subroutine duscr_begin(m, n, A, istat)
Allocates an empty matrix (A) and leaves it in build state.
Definition: rsb_blas_sparse.F90:375
Definition: rsb.F90:11
integer, parameter rsb_idx_kind
Definition: rsb.F90:15
integer(c_signed_char), parameter rsb_numerical_type_double_complex
Definition: rsb.F90:1081
integer(c_int), parameter rsb_io_want_extra_verbose_interface
See RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE.
Definition: rsb.F90:1111
integer(c_int), parameter rsb_io_want_verbose_tuning
See RSB_IO_WANT_VERBOSE_TUNING.
Definition: rsb.F90:1138

Function/Subroutine Documentation

◆ blas_sparse_io_example()

subroutine blas_sparse_io_example ( integer(kind=rsb_idx_kind)  res)

◆ blas_sparse_mod_example()

subroutine blas_sparse_mod_example ( integer(kind=rsb_blas_idx_kind)  res)

◆ main()

program main ( void  )