-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMODULE BASIS.f90
318 lines (298 loc) · 14.7 KB
/
MODULE BASIS.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
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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
!***********************************************************
!* Fortran 90 version of basis_r.cpp, reduced and modified *
!* version of basis.c by J-P Moreau, Paris *
!* ------------------------------------------------------- *
!* Reference for original basis.c: *
!* *
!* "Numerical Algorithms with C, by Gisela Engeln-Muellges *
!* and Frank Uhlig, Springer-Verlag, 1996". *
!***********************************************************
MODULE BASIS
PARAMETER(HALF=0.5d0,ONE=1.d0,TWO=2.d0,THREE=3.d0)
CONTAINS
Subroutine NormMax(vektor, n, normaxi)
DOUBLE PRECISION normaxi, vektor(0:n-1)
!************************************************************************
!* Return the maximum norm of a (0..n-1) vector v. *
!* *
!* global names used: *
!* ================== *
!* None. *
!************************************************************************
DOUBLE PRECISION norm ! local max
DOUBLE PRECISION betrag ! magnitude of a component
norm=0.d0
do i=0, n-1
betrag=abs(vektor(i))
if (betrag > norm) norm = betrag
end do
normaxi = norm
return
end subroutine NormMax
Subroutine CopyMat(ziel, quelle, n, m)
DOUBLE PRECISION ziel(0:n-1,0:m-1),quelle(0:n-1,0:m-1)
!************************************************************************
!* copy the n*m elements of the matrix quelle into the matrix ziel. *
!* *
!* global names used: *
!* ================== *
!* none *
!************************************************************************
do i=0, n-1
do j=0, m-1
ziel(i,j) = quelle(i,j)
end do
end do
return
end subroutine CopyMat
!*****************************************
! Read a vector from unit u from index = 0
!*****************************************
Subroutine ReadVec (u, n, x)
INTEGER u
DOUBLE PRECISION x(0:n-1)
read(u,*) (x(i),i=0,n-1)
return
end subroutine ReadVec
!*****************************************
! Read a vector from unit u from index = 1
!*****************************************
Subroutine ReadVec1 (u, n, x)
INTEGER u
DOUBLE PRECISION x(1:n)
read(u,*) (x(i),i=1,n)
return
end subroutine ReadVec1
!*****************************
! Put all the components of a
! vector to a given value val
!*****************************
Subroutine SetVec(n,x,val)
DOUBLE PRECISION x(0:n-1)
DOUBLE PRECISION val
do i = 0, n-1
x(i) = val
end do
return
end subroutine SetVec
!********************************
! Put all the components of a
! nxm matrix to a given value val
!********************************
Subroutine SetMat(n,m,a,val)
DOUBLE PRECISION a(0:n-1,0:m-1)
DOUBLE PRECISION val
do i = 0, n-1
do j = 0, m-1
a(i,j) = val
end do
end do
return
end subroutine SetMat
Subroutine WriteVec(u, n, x)
INTEGER u, n
DOUBLE PRECISION x(0:n-1)
!*====================================================================*
!* *
!* Put out vector x of length n to output file (10 items per line) *
!* Index starts at ZERO. *
!* *
!*====================================================================*
!* *
!* Input parameters: *
!* ================ *
!* *
!* n integer n; *
!* lenght of vector *
!* x DOUBLE PRECISION*8 x(n) *
!* vector *
!* *
!*====================================================================*
write(u,10) (x(i),i=0,n-1)
return
10 format(10E13.6)
end subroutine WriteVec
Subroutine WriteVec1 (u, n, x)
INTEGER u, n
DOUBLE PRECISION x(1:n)
!*====================================================================*
!* *
!* Put out vector x of length n to output file (10 items per line) *
!* Index starts at ONE. *
!* *
!*====================================================================*
!* *
!* Input parameters: *
!* ================ *
!* *
!* n integer n; *
!* lenght of vector *
!* x DOUBLE PRECISION*8 x(n) *
!* vector *
!* *
!*====================================================================*
write(u,10) (x(i),i=1,n)
return
10 format(10E13.6)
end subroutine WriteVec1
Subroutine ReadMat (u,n, m, a)
!*====================================================================*
!* *
!* Read an n x m matrix from input file. *
!* *
!*====================================================================*
!* *
!* Input parameters : *
!* ================== *
!* u integer u; # of unit to read *
!* n integer m; ( m > 0 ) *
!* number of rows of matrix *
!* m integer n; ( n > 0 ) *
!* column number of matrix *
!* *
!* Output parameter: *
!* ================ *
!* a DOUBLE PRECISION*8 a(n,m) *
!* matrix to read *
!* *
!* ATTENTION : WE do not allocate storage for a here. *
!* *
!*====================================================================*
INTEGER u, n, m
DOUBLE PRECISION a(0:n-1,0:m-1)
Read(u,*) ((a(i,j),j=0,m-1),i=0,n-1)
return
end subroutine ReadMat
Subroutine ReadMat1(u,n, m, a)
!*====================================================================*
!* *
!* Read an n x m matrix from input file (index begins at 1) *
!* *
!*====================================================================*
!* *
!* Input parameters : *
!* ================== *
!* u integer u; # of unit to read *
!* n integer m; ( m > 0 ) *
!* number of rows of matrix *
!* m integer n; ( n > 0 ) *
!* column number of matrix *
!* *
!* Output parameter: *
!* ================ *
!* a DOUBLE PRECISION*8 a(n,m) *
!* matrix to read *
!* *
!* ATTENTION : WE do not allocate storage for a here. *
!* *
!*====================================================================*
INTEGER u, n, m
DOUBLE PRECISION a(n,m)
Read(u,*) ((a(i,j),j=1,m),i=1,n)
return
end subroutine ReadMat1
Subroutine WriteMat(u,n, m, a)
!*====================================================================*
!* *
!* Put out an m x n matrix in output file ( 10 items per line ) *
!* *
!*====================================================================*
!* *
!* Input parameters: *
!* ================ *
!* u integer u; # of unit to read *
!* n int m; ( m > 0 ) *
!* row number of matrix *
!* m int n; ( n > 0 ) *
!* column number of matrix *
!* a DOUBLE PRECISION*8 a(n,n) *
!* matrix to write *
!*====================================================================*
INTEGER u
DOUBLE PRECISION a(0:n-1,0:m-1)
do i=0, n-1
write(u,10) (a(i,j),j=0,m-1)
enddo
return
10 format(10E14.6)
end subroutine WriteMat
Subroutine WriteMat1(u,n, m, a)
!*====================================================================*
!* *
!* Put out an m x n matrix in output file ( 10 items per line ) *
!* (index begins at one) *
!* *
!*====================================================================*
!* *
!* Input parameters: *
!* ================ *
!* u integer u; # of unit to read *
!* n int m; ( m > 0 ) *
!* row number of matrix *
!* m int n; ( n > 0 ) *
!* column number of matrix *
!* a DOUBLE PRECISION*8 a(n,n) *
!* matrix to write *
!*====================================================================*
INTEGER u
DOUBLE PRECISION a(n,m)
do i=1, n
write(u,10) (a(i,j),j=1,m)
enddo
return
10 format(10E14.6)
end subroutine WriteMat1
Subroutine WriteHead (u, nom)
!*====================================================================*
!* *
!* Put out header with text in unit u . *
!* *
!*====================================================================*
!* *
!* Input parameters: *
!* ================ *
!* u integer u; # of unit to read *
!* nom character*(*) nom *
!* text of headertext (last byte is a 0) *
!* *
!* Return value : *
!* ============= *
!* None *
!* *
!*====================================================================*
INTEGER u
CHARACTER*(*) nom
Write(u,*) '----------------------------------------------------------'
Write(u,*) nom
Write(u,*) '----------------------------------------------------------'
return
end subroutine WriteHead
Subroutine WriteEnd(u)
!*====================================================================*
!* *
!* Put out end of writing onto unit u. *
!* *
!*====================================================================*
!* *
!* Return value : *
!* ============= *
!* None *
!* *
!*====================================================================*
INTEGER u
Write(u,*) ' '
Write(u,*) '----------------------------------------------------------'
return
end subroutine WriteEnd
!***********************
! Swap two DOUBLE PRECISION values
!***********************
Subroutine Swap(a,b)
DOUBLE PRECISION a,b,temp
temp=a
a=b
b=temp
return
end subroutine Swap
END MODULE
! end of file basis.f90