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
|
! A simple OpenMP example to use SuperLU to solve multiple independent linear systems.
! Contributor: Ed D'Azevedo, Oak Ridge National Laboratory
!
program tslu_omp
implicit none
integer, parameter :: maxn = 10*1000
integer, parameter :: maxnz = 100*maxn
integer, parameter :: nsys = 6 !! 64
real*8 :: values(maxnz), b(maxn)
integer :: rowind(maxnz), colptr(maxn)
! integer :: Ai(maxnz, nsys), Aj(maxn, nsys) ! Sherry added
integer n, nnz, nrhs, ldb, info, iopt
integer*8 :: factors, lufactors(nsys)
real*8 :: A(maxnz, nsys)
integer :: luinfo(nsys)
real*8 :: brhs(maxn,nsys)
integer :: i,j
real*8 :: err, maxerr
integer :: nthread
!$ integer, external :: omp_get_num_threads
! --------------
! read in matrix
! --------------
print*, 'before hbcode1'
call hbcode1(n,n,nnz,values,rowind,colptr)
print*, 'after hbcode1'
nthread = 1
!$omp parallel
!$omp master
!$ nthread = omp_get_num_threads()
!$omp end master
!$omp end parallel
write(*,*) 'nthreads = ',nthread
write(*,*) 'nsys = ',nsys
write(*,*) 'n, nnz ', n, nnz
!$omp parallel do private(j)
do j=1,nsys
A(1:nnz,j) = values(1:nnz)
enddo
nrhs = 1
ldb = n
!$omp parallel do private(j)
do j=1,nsys
brhs(:,j) = j
enddo
! ---------------------
! perform factorization
! ---------------------
iopt = 1
!$omp parallel do private(j,values,b,info,factors)
do j=1,nsys
!$omp parallel workshare
values(1:nnz) = A(1:nnz,j)
b(1:n) = brhs(1:n,j)
!$omp end parallel workshare
info = 0
call c_fortran_dgssv( iopt,n,nnz, nrhs, values, rowind, colptr, &
& b, ldb, factors, info )
!$omp parallel workshare
A(1:nnz,j) = values(1:nnz)
brhs(1:n,j) = b(1:n)
!$omp end parallel workshare
luinfo(j) = info
lufactors(j) = factors
enddo
do j=1,nsys
info = luinfo(j)
if (info.ne.0) then
write(*,9010) j, info
9010 format(' factorization of j=',i7,' returns info= ',i7)
endif
enddo
! ---------------------------------------
! solve the system using existing factors
! ---------------------------------------
iopt = 2
!$omp parallel do private(j,b,values,factors,info)
do j=1,nsys
factors = lufactors(j)
values(1:nnz) = A(1:nnz,j)
info = 0
b(1:n) = brhs(1:n,j)
call c_fortran_dgssv( iopt,n,nnz,nrhs,values,rowind,colptr, &
& b,ldb,factors,info )
lufactors(j) = factors
luinfo(j) = info
brhs(1:n,j) = b(1:n)
enddo
! ------------
! simple check
! ------------
err = 0
maxerr = 0
do j=2,nsys
do i=1,n
err = abs(brhs(i,1)*j - brhs(i,j))
maxerr = max(maxerr,err)
enddo
enddo
write(*,*) 'max error = ', maxerr
! -------------
! free storage
! -------------
iopt = 3
!$omp parallel do private(j)
do j=1,nsys
call c_fortran_dgssv(iopt,n,nnz,nrhs,A(:,j),rowind,colptr, &
& brhs(:,j), ldb, lufactors(j), luinfo(j) )
enddo
stop
end program
|
