-
Notifications
You must be signed in to change notification settings - Fork 1
/
setuprhs.f90
88 lines (54 loc) · 1.33 KB
/
setuprhs.f90
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
subroutine setupRHS(stat,dim,npts,fpcb,gpcb,hpcb,numpc,rhsF)
implicit none
include 'collsub.h'
integer,intent(in):: stat,dim,npts
integer::j,k,i
integer,intent(out)::numpc
real*8,intent(in) :: fpcb(MAXPTS),gpcb(DIM,MAXPTS),hpcb(dim,dim,maxpts)
real*8::rhsF(maxdat)
real*8 :: tmp(maxdat)
if (stat.eq.0) then
numpc=0
do j=1,npts
numpc=numpc+1
tmp(j) = fpcb(j)
enddo
rhsF(:)=tmp(:)
else if (stat.eq.1) then
numpc=0
do j=1,npts
numpc=numpc+1
tmp(numpc) = fpcb(j)
do k=1,DIM
numpc=numpc+1
tmp(numpc) = gpcb(k,j)
end do
enddo
rhsF=tmp
else if (stat.eq.2) then
numpc=0
do j=1,npts
numpc=numpc+1
tmp(numpc) = fpcb(j)
do k=1,DIM
numpc=numpc+1
tmp(numpc) = gpcb(k,j)
end do
do k=1,dim
do i=1,dim
if (k.eq.i.or.k.gt.i) numpc=numpc+1
if (k.eq.i) then !diagonal
tmp(numpc) = hpcb(k,i,j)
else if (k.gt.i) then !upper diagonal
tmp(numpc) = hpcb(k,i,j)
else
! numpc=numpc-1
end if
end do
end do
enddo
rhsF=tmp
else
stop 'unknown stat'
end if
end subroutine setupRHS