1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
|
!
! Copyright (C) 2008-2016 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/>.
!
SUBROUTINE rsb_mod_example1(res)
USE rsb
USE ISO_C_BINDING
IMPLICIT NONE
INTEGER ::res
INTEGER,TARGET :: istat = 0, i
INTEGER :: transt = RSB_TRANSPOSITION_N ! Please note that this interface is unfinished
INTEGER :: incx = 1, incy = 1
REAL(KIND=8),TARGET :: alpha = 3, beta = 1
! 1 1
! 1 1
! declaration of VA,IA,JA
INTEGER :: nnz = 4
INTEGER :: nr = 2
INTEGER :: nc = 2
INTEGER :: nrhs = 1
INTEGER :: order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER ! rhs layout
INTEGER :: flags = RSB_FLAG_NOFLAGS
INTEGER,TARGET :: IA(4) = (/0, 1, 1,0/)
INTEGER,TARGET :: JA(4) = (/0, 0, 1,1/)
REAL(KIND=8),TARGET :: VA(4) = (/1,1,1,1/)
REAL(KIND=8),TARGET :: x(2) = (/1, 1/)! reference x
REAL(KIND=8),TARGET :: cy(2) = (/9, 9/)! reference cy after
REAL(KIND=8),TARGET :: y(2) = (/3, 3/)! y will be overwritten
TYPE(C_PTR),TARGET :: mtxAp = C_NULL_PTR ! matrix pointer
REAL(KIND=8) :: tmax = 2.0 ! tuning max time
INTEGER :: titmax = 2 ! tuning max iterations
INTEGER,TARGET :: ont = 0 ! optimal number of threads
res = 0
mtxAp = rsb_mtx_alloc_from_coo_const(C_LOC(VA),C_LOC(IA),C_LOC(JA)&
&,nnz,&
& RSB_NUMERICAL_TYPE_DOUBLE,nr,nc,1,1,flags,C_LOC(istat))
IF (istat.NE.RSB_ERR_NO_ERROR) GOTO 9997
istat = rsb_file_mtx_save(mtxAp,C_NULL_PTR)
! Structure autotuning:
istat = rsb_tune_spmm(C_LOC(mtxAp),C_NULL_PTR,C_NULL_PTR,titmax,&
& tmax,&
& transt,C_LOC(alpha),C_NULL_PTR,nrhs,order,C_LOC(x),nr,&
& C_LOC(beta),C_LOC(y),nc)
IF (istat.NE.RSB_ERR_NO_ERROR) GOTO 9997
! Thread count autotuning:
istat = rsb_tune_spmm(C_NULL_PTR,C_NULL_PTR,C_LOC(ont),titmax,&
& tmax,&
& transt,C_LOC(alpha),mtxAp,nrhs,order,C_LOC(x),nr,C_LOC(beta),&
& C_LOC(y),nc)
PRINT *, "Optimal number of threads:", ont
y(:) = (/3, 3/)! reference y
IF (istat.NE.RSB_ERR_NO_ERROR) GOTO 9997
istat = rsb_file_mtx_save(mtxAp,C_NULL_PTR)
IF (istat.NE.RSB_ERR_NO_ERROR) GOTO 9997
istat = rsb_spmv(transt,C_LOC(alpha),mtxAp,C_LOC(x),incx,&
& C_LOC(beta),C_LOC(y),incy)
IF (istat.NE.RSB_ERR_NO_ERROR) GOTO 9997
DO i = 1, 2
IF (y(i).NE.cy(i)) PRINT *, "type=d dims=2x2 sym=g diag=g &
&blocks=1x1 usmv alpha= 3 beta= 1 incx=1 incy=1 trans=n is not ok"
IF (y(i).NE.cy(i)) GOTO 9997
END DO
PRINT*,"type=d dims=2x2 sym=g diag=g blocks=1x1 usmv alpha= 3&
& beta= 1 incx=1 incy=1 trans=n is ok"
GOTO 9998
9997 res = -1
9998 CONTINUE
mtxAp = rsb_mtx_free(mtxAp)
IF (istat.NE.RSB_ERR_NO_ERROR) res = -1
! 9999 CONTINUE
istat = rsb_perror(C_NULL_PTR,istat)
end SUBROUTINE rsb_mod_example1
SUBROUTINE rsb_mod_example2(res)
USE rsb
USE ISO_C_BINDING
IMPLICIT NONE
INTEGER,TARGET :: errval
INTEGER :: res
INTEGER :: transt = RSB_TRANSPOSITION_N ! no transposition
INTEGER :: incX = 1, incB = 1 ! X, B vectors increment
REAL(KIND=8),TARGET :: alpha = 3,beta = 1
INTEGER :: nnzA = 4, nrA = 3, ncA = 3 ! nonzeroes, rows, columns of matrix A
INTEGER,TARGET :: IA(4) = (/1, 2, 3, 3/) ! row indices
INTEGER,TARGET :: JA(4) = (/1, 2, 1, 3/) ! column indices
INTEGER(C_SIGNED_CHAR) :: typecode = RSB_NUMERICAL_TYPE_DOUBLE
INTEGER :: flags =RSB_FLAG_DEFAULT_MATRIX_FLAGS+RSB_FLAG_SYMMETRIC
REAL(KIND=8),TARGET :: VA(4) = (/11.0, 22.0, 13.0, 33.0/) ! coefficients
REAL(KIND=8),TARGET :: X(3) = (/ 0, 0, 0/)
REAL(KIND=8),TARGET :: B(3) = (/-1.0, -2.0, -2.0/)
TYPE(C_PTR),TARGET :: mtxAp = C_NULL_PTR
TYPE(C_PTR) :: mtxApp = C_NULL_PTR
REAL(KIND=8),TARGET :: ETIME = 0.0
!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
errval = rsb_lib_init(IO) ! librsb initialization
IF (errval.NE.RSB_ERR_NO_ERROR) &
& STOP "error calling rsb_lib_init"
#if defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ < 5)
#define RSB_SKIP_BECAUSE_OLD_COMPILER 1
#endif
#ifndef RSB_SKIP_BECAUSE_OLD_COMPILER
mtxAp = rsb_mtx_alloc_from_coo_begin(nnzA,typecode,nrA,ncA,flags,&
& C_LOC(errval)) ! begin matrix creation
errval = rsb_mtx_set_vals(mtxAp,&
& C_LOC(VA),C_LOC(IA),C_LOC(JA),nnzA,flags) ! insert some nonzeroes
mtxApp = C_LOC(mtxAp) ! Old compilers like e.g.: Gfortran 4.4.7 will NOT compile this.
IF (errval.NE.RSB_ERR_NO_ERROR) &
& STOP "error calling rsb_mtx_set_vals"
errval = rsb_mtx_alloc_from_coo_end(mtxApp) ! end matrix creation
IF (errval.NE.RSB_ERR_NO_ERROR) &
& STOP "error calling rsb_mtx_alloc_from_coo_end"
errval = rsb_spmv(transt,C_LOC(alpha),mtxAp,C_LOC(X),&
& incX,C_LOC(beta),C_LOC(B),incB) ! X := X + (3) * A * B
IF (errval.NE.RSB_ERR_NO_ERROR)&
& STOP "error calling rsb_spmv"
mtxAp = rsb_mtx_free(mtxAp) ! destroy matrix
! The following is optional and depends on configure options, so it is allowed to fail
errval = rsb_lib_get_opt(RSB_IO_WANT_LIBRSB_ETIME,C_LOC(ETIME))
IF (errval.EQ.RSB_ERR_NO_ERROR)&
& PRINT*,"Time spent in librsb is:",ETIME
! IF (errval.NE.0)STOP "error calling rsb_lib_get_opt"
errval = RSB_ERR_NO_ERROR
IF (errval.NE.RSB_ERR_NO_ERROR) &
& STOP "error calling rsb_mtx_free"
#else
PRINT*,"You have an old Fortran compiler not supporting C_LOC."
PRINT*,"Skipping a part of the test"
#endif
errval=rsb_lib_exit(EO) ! librsb finalization
IF (errval.NE.RSB_ERR_NO_ERROR)&
& STOP "error calling rsb_lib_exit"
PRINT *, "rsb module fortran test is ok"
res = errval
end SUBROUTINE rsb_mod_example2
PROGRAM main
USE rsb
IMPLICIT NONE
INTEGER :: res = RSB_ERR_NO_ERROR, passed = 0, failed = 0
!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
res = rsb_lib_init(IO)
CALL rsb_mod_example1(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
res = rsb_lib_exit(EO)
CALL rsb_mod_example2(res)
IF (res.LT.0) failed = failed + 1
IF (res.EQ.0) passed = passed + 1
PRINT *, "FAILED:", failed
PRINT *, "PASSED:", passed
IF (failed.GT.0) THEN
STOP 1
END IF
END PROGRAM
|
