CoCalc -- huti

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fhutiter/src/huti_aux.F90
³²⁰⁶ views
1
module huti_aux
2

3
  !
4
  ! *
5
  ! *  Elmer, A Finite Element Software for Multiphysical Problems
6
  ! *
7
  ! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
8
  ! * 
9
  ! * This library is free software; you can redistribute it and/or
10
  ! * modify it under the terms of the GNU Lesser General Public
11
  ! * License as published by the Free Software Foundation; either
12
  ! * version 2.1 of the License, or (at your option) any later version.
13
  ! *
14
  ! * This library is distributed in the hope that it will be useful,
15
  ! * but WITHOUT ANY WARRANTY; without even the implied warranty of
16
  ! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
  ! * Lesser General Public License for more details.
18
  ! * 
19
  ! * You should have received a copy of the GNU Lesser General Public
20
  ! * License along with this library (in file ../LGPL-2.1); if not, write 
21
  ! * to the Free Software Foundation, Inc., 51 Franklin Street, 
22
  ! * Fifth Floor, Boston, MA  02110-1301  USA
23
  ! *
24
  ! **************************************************************************/
25

26
  ! Auxiliary routines for HUTI
27
  !
28

29
#include  "huti_fdefs.h" 
30

31
contains
32

33
  !*************************************************************************
34
  !*************************************************************************
35
  !
36
  ! This is a dummy preconditioning routine copying only the vector
37

38
  subroutine  huti_sdummy_pcondfun  ( u, v, ipar )
39

40
    implicit none
41

42
    real, dimension(*) :: u, v
43
    integer, dimension(*) :: ipar
44

45
    integer :: i
46

47
    !*************************************************************************
48

49
    do i = 1, HUTI_NDIM
50
       u(i) = v(i)
51
    end do
52

53
    return
54

55
  end subroutine  huti_sdummy_pcondfun
56

57
  !*************************************************************************
58

59
  !*************************************************************************
60
  !*************************************************************************
61
  !
62
  ! This routine fills a vector with pseudo random numbers
63

64
  subroutine  huti_srandvec  ( u, ipar )
65

66
    implicit none
67

68
    real, dimension(*) :: u
69
    integer, dimension(*) :: ipar
70

71
    integer :: i
72
    real :: harvest
73

74
    !*************************************************************************
75

76
    do i = 1, HUTI_NDIM
77
       call random_number( harvest )
78
       u(i) = harvest
79
    end do
80

81
    return
82

83
  end subroutine  huti_srandvec
84

85
  !*************************************************************************
86
  !*************************************************************************
87
  !
88
  ! Construct LU decomposition of the given matrix and solve LUu = v
89
  !
90

91
  subroutine  huti_slusolve  ( n, lumat, u, v )
92

93
    implicit none
94

95
    integer :: n
96
    real, dimension(n,n) :: lumat
97
    real, dimension(n) :: u, v
98

99
    integer :: i, j, k
100

101
    !*************************************************************************
102

103
    !
104
    ! This is from Saad''s book, Algorithm 10.4
105
    !
106

107
    do i = 2,n
108
       do k = 1, i - 1
109

110
          ! Check for small pivot
111

112
          if ( abs(lumat(k,k)) .lt. 1.0e-16 ) then
113
             print *, '(libhuti.a) GMRES: small pivot', lumat(k,k)
114
          end if
115

116
          ! Compute a_ik = a_ik / a_kk
117

118
          lumat(i,k) = lumat(i,k) / lumat(k,k)
119

120
          do j = k + 1, n
121

122
             ! Compute a_ij = a_ij - a_ik * a_kj
123

124
             lumat(i,j) = lumat(i,j) - lumat(i,k) * lumat(k,j)
125
          end do
126

127
       end do
128
    end do
129

130
    ! Forward solve, Lu = v
131

132
    do i = 1, n
133

134
       ! Compute u(i) = v(i) - sum L(i,j) u(j)
135

136
       u(i) = v(i)
137
       do k = 1, i - 1
138
          u(i) = u(i) - lumat(i,k) * u(k)
139
       end do
140
    end do
141

142
    ! Backward solve, u = inv(U) u
143

144
    do i = n, 1, -1
145

146
       ! Compute u(i) = u(i) - sum U(i,j) u(j)
147

148
       do k = i + 1, n
149
          u(i) = u(i) - lumat(i,k) * u(k)
150
       end do
151

152
       ! Compute u(i) = u(i) / U(i,i)
153

154
       u(i) = u(i) / lumat(i,i) 
155
    end do
156

157
    return
158

159
  end subroutine  huti_slusolve
160

161

162
  !*************************************************************************
163
  !*************************************************************************
164
  !
165
  ! This is a dummy preconditioning routine copying only the vector
166

167
  subroutine  huti_ddummy_pcondfun  ( u, v, ipar )
168

169
    implicit none
170

171
    double precision, dimension(*) :: u, v
172
    integer, dimension(*) :: ipar
173

174
    integer :: i
175

176
    !*************************************************************************
177
    !$OMP PARALLEL DO
178
    do i = 1, HUTI_NDIM
179
       u(i) = v(i)
180
    end do
181
    !$OMP END PARALLEL DO 
182

183
    return
184

185
  end subroutine  huti_ddummy_pcondfun
186

187
  !*************************************************************************
188

189
  !*************************************************************************
190
  !*************************************************************************
191
  !
192
  ! This routine fills a vector with pseudo random numbers
193

194
  subroutine  huti_drandvec  ( u, ipar )
195

196
    implicit none
197

198
    double precision, dimension(*) :: u
199
    integer, dimension(*) :: ipar
200

201
    integer :: i
202
    real :: harvest
203

204
    !*************************************************************************
205

206
    do i = 1, HUTI_NDIM
207
       call random_number( harvest )
208
       u(i) = harvest
209
    end do
210

211
    return
212

213
  end subroutine  huti_drandvec
214

215
  !*************************************************************************
216
  !*************************************************************************
217
  !
218
  ! Construct LU decomposition of the given matrix and solve LUu = v
219
  !
220

221
  subroutine  huti_dlusolve  ( n, lumat, u, v )
222

223
    implicit none
224

225
    integer :: n
226
    double precision, dimension(n,n) :: lumat
227
    double precision, dimension(n) :: u, v
228

229
    integer :: i, j, k
230

231
    !*************************************************************************
232

233
    !
234
    ! This is from Saad''s book, Algorithm 10.4
235
    !
236

237
    do i = 2,n
238
       do k = 1, i - 1
239

240
          ! Check for small pivot
241

242
          if ( abs(lumat(k,k)) .lt. 1.0e-16 ) then
243
             print *, '(libhuti.a) GMRES: small pivot', lumat(k,k)
244
          end if
245

246
          ! Compute a_ik = a_ik / a_kk
247

248
          lumat(i,k) = lumat(i,k) / lumat(k,k)
249

250
          do j = k + 1, n
251

252
             ! Compute a_ij = a_ij - a_ik * a_kj
253

254
             lumat(i,j) = lumat(i,j) - lumat(i,k) * lumat(k,j)
255
          end do
256

257
       end do
258
    end do
259

260
    ! Forward solve, Lu = v
261

262
    do i = 1, n
263

264
       ! Compute u(i) = v(i) - sum L(i,j) u(j)
265

266
       u(i) = v(i)
267
       do k = 1, i - 1
268
          u(i) = u(i) - lumat(i,k) * u(k)
269
       end do
270
    end do
271

272
    ! Backward solve, u = inv(U) u
273

274
    do i = n, 1, -1
275

276
       ! Compute u(i) = u(i) - sum U(i,j) u(j)
277

278
       do k = i + 1, n
279
          u(i) = u(i) - lumat(i,k) * u(k)
280
       end do
281

282
       ! Compute u(i) = u(i) / U(i,i)
283

284
       u(i) = u(i) / lumat(i,i) 
285
    end do
286

287
    return
288

289
  end subroutine  huti_dlusolve
290

291

292
  !*************************************************************************
293
  !*************************************************************************
294
  !
295
  ! This is a dummy preconditioning routine copying only the vector
296

297
  subroutine  huti_cdummy_pcondfun  ( u, v, ipar )
298

299
    implicit none
300

301
    complex, dimension(*) :: u, v
302
    integer, dimension(*) :: ipar
303

304
    integer :: i
305

306
    !*************************************************************************
307

308
    do i = 1, HUTI_NDIM
309
       u(i) = v(i)
310
    end do
311

312
    return
313

314
  end subroutine  huti_cdummy_pcondfun
315

316
  !*************************************************************************
317

318
  !*************************************************************************
319
  !*************************************************************************
320
  !
321
  ! This routine fills a vector with pseudo random numbers
322

323
  subroutine  huti_crandvec  ( u, ipar )
324

325
    implicit none
326

327
    complex, dimension(*) :: u
328
    integer, dimension(*) :: ipar
329

330
    integer :: i
331
    real :: harvest
332

333
    !*************************************************************************
334

335
    do i = 1, HUTI_NDIM
336
       call random_number( harvest )
337
       u(i) = harvest
338
       call random_number( harvest )
339
       u(i) = cmplx( 0, harvest )
340
    end do
341

342
    return
343

344
  end subroutine  huti_crandvec
345

346
  !*************************************************************************
347
  !*************************************************************************
348
  !
349
  ! Construct LU decomposition of the given matrix and solve LUu = v
350
  !
351

352
  subroutine  huti_clusolve  ( n, lumat, u, v )
353

354
    implicit none
355

356
    integer :: n
357
    complex, dimension(n,n) :: lumat
358
    complex, dimension(n) :: u, v
359

360
    integer :: i, j, k
361

362
    !*************************************************************************
363

364
    !
365
    ! This is from Saad''s book, Algorithm 10.4
366
    !
367

368
    do i = 2,n
369
       do k = 1, i - 1
370

371
          ! Check for small pivot
372

373
          if ( abs(lumat(k,k)) .lt. 1.0e-16 ) then
374
             print *, '(libhuti.a) GMRES: small pivot', lumat(k,k)
375
          end if
376

377
          ! Compute a_ik = a_ik / a_kk
378

379
          lumat(i,k) = lumat(i,k) / lumat(k,k)
380

381
          do j = k + 1, n
382

383
             ! Compute a_ij = a_ij - a_ik * a_kj
384

385
             lumat(i,j) = lumat(i,j) - lumat(i,k) * lumat(k,j)
386
          end do
387

388
       end do
389
    end do
390

391
    ! Forward solve, Lu = v
392

393
    do i = 1, n
394

395
       ! Compute u(i) = v(i) - sum L(i,j) u(j)
396

397
       u(i) = v(i)
398
       do k = 1, i - 1
399
          u(i) = u(i) - lumat(i,k) * u(k)
400
       end do
401
    end do
402

403
    ! Backward solve, u = inv(U) u
404

405
    do i = n, 1, -1
406

407
       ! Compute u(i) = u(i) - sum U(i,j) u(j)
408

409
       do k = i + 1, n
410
          u(i) = u(i) - lumat(i,k) * u(k)
411
       end do
412

413
       ! Compute u(i) = u(i) / U(i,i)
414

415
       u(i) = u(i) / lumat(i,i) 
416
    end do
417

418
    return
419

420
  end subroutine  huti_clusolve
421

422

423
  !*************************************************************************
424
  !*************************************************************************
425
  !
426
  ! This is a dummy preconditioning routine copying only the vector
427

428
  subroutine  huti_zdummy_pcondfun  ( u, v, ipar )
429

430
    implicit none
431

432
    double complex, dimension(*) :: u, v
433
    integer, dimension(*) :: ipar
434

435
    integer :: i
436

437
    !*************************************************************************
438

439
    do i = 1, HUTI_NDIM
440
       u(i) = v(i)
441
    end do
442

443
    return
444

445
  end subroutine  huti_zdummy_pcondfun
446

447
  !*************************************************************************
448

449
  !*************************************************************************
450
  !*************************************************************************
451
  !
452
  ! This routine fills a vector with pseudo random numbers
453

454
  subroutine  huti_zrandvec  ( u, ipar )
455

456
    implicit none
457

458
    double complex, dimension(*) :: u
459
    integer, dimension(*) :: ipar
460

461
    integer :: i
462
    real :: harvest
463

464
    !*************************************************************************
465

466
    do i = 1, HUTI_NDIM
467
       call random_number( harvest )
468
       u(i) = harvest
469
       call random_number( harvest )
470
       u(i) = cmplx( 0, harvest )
471
    end do
472

473
    return
474

475
  end subroutine  huti_zrandvec
476

477
  !*************************************************************************
478
  !*************************************************************************
479
  !
480
  ! Construct LU decomposition of the given matrix and solve LUu = v
481
  !
482

483
  subroutine  huti_zlusolve  ( n, lumat, u, v )
484

485
    implicit none
486

487
    integer :: n
488
    double complex, dimension(n,n) :: lumat
489
    double complex, dimension(n) :: u, v
490

491
    integer :: i, j, k
492

493
    !*************************************************************************
494

495
    !
496
    ! This is from Saad''s book, Algorithm 10.4
497
    !
498

499
    do i = 2,n
500
       do k = 1, i - 1
501

502
          ! Check for small pivot
503

504
          if ( abs(lumat(k,k)) .lt. 1.0e-16 ) then
505
             print *, '(libhuti.a) GMRES: small pivot', lumat(k,k)
506
          end if
507

508
          ! Compute a_ik = a_ik / a_kk
509

510
          lumat(i,k) = lumat(i,k) / lumat(k,k)
511

512
          do j = k + 1, n
513

514
             ! Compute a_ij = a_ij - a_ik * a_kj
515

516
             lumat(i,j) = lumat(i,j) - lumat(i,k) * lumat(k,j)
517
          end do
518

519
       end do
520
    end do
521

522
    ! Forward solve, Lu = v
523

524
    do i = 1, n
525

526
       ! Compute u(i) = v(i) - sum L(i,j) u(j)
527

528
       u(i) = v(i)
529
       do k = 1, i - 1
530
          u(i) = u(i) - lumat(i,k) * u(k)
531
       end do
532
    end do
533

534
    ! Backward solve, u = inv(U) u
535

536
    do i = n, 1, -1
537

538
       ! Compute u(i) = u(i) - sum U(i,j) u(j)
539

540
       do k = i + 1, n
541
          u(i) = u(i) - lumat(i,k) * u(k)
542
       end do
543

544
       ! Compute u(i) = u(i) / U(i,i)
545

546
       u(i) = u(i) / lumat(i,i) 
547
    end do
548

549
    return
550

551
  end subroutine  huti_zlusolve
552

553
end module huti_aux
554

555
Product

Resources

Company
