CoCalc -- DirectSolve.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/DirectSolve.F90
⁵²⁰⁹ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer, A Finite Element Software for Multiphysical Problems
4
! *
5
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6
! * 
7
! *  This library is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU Lesser General Public
9
! *  License as published by the Free Software Foundation; either
10
! *  version 2.1 of the License, or (at your option) any later version.
11
! *
12
! *  This library is distributed in the hope that it will be useful,
13
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
! *  Lesser General Public License for more details.
16
! * 
17
! *  You should have received a copy of the GNU Lesser General Public
18
! *  License along with this library (in file ../LGPL-2.1); if not, write 
19
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
20
! *  Fifth Floor, Boston, MA  02110-1301  USA
21
! *
22
! *****************************************************************************/
23
!
24
!/******************************************************************************
25
! *
26
! *  Authors: Juha Ruokolainen
27
! *  Email:   [email protected]
28
! *  Web:     http://www.csc.fi/elmer
29
! *  Address: CSC - IT Center for Science Ltd.
30
! *           Keilaranta 14
31
! *           02101 Espoo, Finland 
32
! *
33
! *  Original Date: 08 Jun 1997
34
! *
35
! *****************************************************************************/
36

37
!> \ingroup ElmerLib 
38
!> \{
39

40
!------------------------------------------------------------------------------
41
!>  Module containing the direct solvers for linear systems given in CRS format.
42
!> Included are Lapack band matrix solver, multifrontal Umfpack, MUMPS, SuperLU, 
43
!> and Pardiso. Note that many of these are linked in with ElmerSolver only 
44
!> if they are made available at the compilation time. 
45
!------------------------------------------------------------------------------
46

47
#include "../config.h"
48

49
MODULE DirectSolve
50

51
   USE CRSMatrix
52
   USE BandMatrix
53
   USE SParIterSolve
54

55
   IMPLICIT NONE
56

57
CONTAINS
58

59

60
!------------------------------------------------------------------------------
61
!> Solver the complex linear system using direct band matrix solver from of Lapack.
62
!------------------------------------------------------------------------------
63
   SUBROUTINE ComplexBandSolver( A,x,b, Free_fact )
64
!------------------------------------------------------------------------------
65

66
     LOGICAL, OPTIONAL :: Free_Fact
67
     TYPE(Matrix_t) :: A
68
     REAL(KIND=dp) :: x(*),b(*)
69
!------------------------------------------------------------------------------
70

71
   
72
     INTEGER :: i,j,k,istat,Subband,N
73
     COMPLEX(KIND=dp), ALLOCATABLE :: BA(:,:)
74

75
     REAL(KIND=dp), POINTER CONTIG :: Values(:)
76
     INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:)
77

78
     SAVE BA
79
!------------------------------------------------------------------------------
80

81
     IF ( PRESENT(Free_Fact) ) THEN
82
       IF ( Free_Fact ) THEN
83
         IF ( ALLOCATED(BA) ) DEALLOCATE(BA)
84
         RETURN
85
       END IF
86
     END IF
87

88
     Rows => A % Rows
89
     Cols => A % Cols
90
     Diag => A % Diag
91
     Values => A % Values
92

93
     n = A % NumberOfRows
94
     x(1:n) = b(1:n)
95
     n = n / 2
96

97
     IF ( A % Format == MATRIX_CRS .AND. .NOT. A % Symmetric ) THEN
98
       Subband = 0
99
       DO i=1,N
100
         DO j=Rows(2*i-1),Rows(2*i)-1,2
101
           Subband = MAX(Subband,ABS((Cols(j)+1)/2-i))
102
         END DO
103
       END DO
104

105
       IF ( .NOT.ALLOCATED( BA ) ) THEN
106

107
         ALLOCATE( BA(3*SubBand+1,N),stat=istat )
108

109
         IF ( istat /= 0 ) THEN
110
           CALL Fatal( 'ComplexBandSolver', 'Memory allocation error.' )
111
         END IF
112

113
       ELSE IF ( SIZE(BA,1) /= 3*Subband+1 .OR. SIZE(BA,2) /= N ) THEN
114

115
         DEALLOCATE( BA )
116
         ALLOCATE( BA(3*SubBand+1,N),stat=istat )
117

118
         IF ( istat /= 0 ) THEN
119
           CALL Fatal( 'ComplexBandSolver', 'Memory allocation error.' )
120
         END IF
121

122
       END IF
123

124
       BA = 0.0D0
125
       DO i=1,N
126
         DO j=Rows(2*i-1),Rows(2*i)-1,2
127
           k = i - (Cols(j)+1)/2 + 2*Subband + 1
128
           BA(k,(Cols(j)+1)/2) = CMPLX(Values(j), -Values(j+1), KIND=dp )
129
         END DO
130
       END DO
131

132
       CALL SolveComplexBandLapack( N,1,BA,x,Subband,3*Subband+1 )
133

134
     ELSE IF ( A % Format == MATRIX_CRS ) THEN
135

136
       Subband = 0
137
       DO i=1,N
138
         DO j=Rows(2*i-1),Diag(2*i-1)
139
           Subband = MAX(Subband,ABS((Cols(j)+1)/2-i))
140
         END DO
141
       END DO
142

143
       IF ( .NOT.ALLOCATED( BA ) ) THEN
144

145
         ALLOCATE( BA(SubBand+1,N),stat=istat )
146

147
         IF ( istat /= 0 ) THEN
148
           CALL Fatal( 'ComplexBandSolver', 'Memory allocation error.' )
149
         END IF
150

151
       ELSE IF ( SIZE(BA,1) /= Subband+1 .OR. SIZE(BA,2) /= N ) THEN
152

153
         DEALLOCATE( BA )
154
         ALLOCATE( BA(SubBand+1,N),stat=istat )
155

156
         IF ( istat /= 0 ) THEN
157
           CALL Fatal( 'ComplexBandSolver', 'Direct solver memory allocation error.' )
158
         END IF
159

160
       END IF
161

162
       BA = 0.0D0
163
       DO i=1,N
164
         DO j=Rows(2*i-1),Diag(2*i-1)
165
           k = i - (Cols(j)+1)/2 + 1
166
           BA(k,(Cols(j)+1)/2) = CMPLX(Values(j), -Values(j+1), KIND=dp )
167
         END DO
168
       END DO
169

170
       CALL SolveComplexSBandLapack( N,1,BA,x,Subband,Subband+1 )
171

172
     END IF
173
!------------------------------------------------------------------------------
174
  END SUBROUTINE ComplexBandSolver 
175
!------------------------------------------------------------------------------
176

177

178
!------------------------------------------------------------------------------
179
!> Solver the real linear system using direct band matrix solver from of Lapack.
180
!------------------------------------------------------------------------------
181
   SUBROUTINE BandSolver( A,x,b,Free_Fact )
182
!------------------------------------------------------------------------------
183
     LOGICAL, OPTIONAL :: Free_Fact
184
     TYPE(Matrix_t) :: A
185
     REAL(KIND=dp) :: x(*),b(*)
186
!------------------------------------------------------------------------------
187

188
     INTEGER :: i,j,k,istat,Subband,N
189
     REAL(KIND=dp), ALLOCATABLE :: BA(:,:)
190

191
     REAL(KIND=dp), POINTER CONTIG :: Values(:)
192
     INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:)
193

194
     SAVE BA
195
!------------------------------------------------------------------------------
196
     IF ( PRESENT(Free_Fact) ) THEN
197
       IF ( Free_Fact ) THEN
198
         IF ( ALLOCATED(BA) ) DEALLOCATE(BA)
199
         RETURN
200
       END IF
201
     END IF
202

203
     N = A % NumberOfRows
204

205
     x(1:n) = b(1:n)
206

207
     Rows => A % Rows
208
     Cols => A % Cols
209
     Diag => A % Diag
210
     Values => A % Values
211

212
     IF ( A % Format == MATRIX_CRS ) THEN ! .AND. .NOT. A % Symmetric ) THEN
213
        Subband = 0
214
        DO i=1,N
215
          DO j=Rows(i),Rows(i+1)-1
216
            Subband = MAX(Subband,ABS(Cols(j)-i))
217
          END DO
218
        END DO
219

220
        IF ( .NOT.ALLOCATED( BA ) ) THEN
221

222
          ALLOCATE( BA(3*SubBand+1,N),stat=istat )
223

224
          IF ( istat /= 0 ) THEN
225
            CALL Fatal( 'BandSolver', 'Memory allocation error.' )
226
          END IF
227

228
        ELSE IF ( SIZE(BA,1) /= 3*Subband+1 .OR. SIZE(BA,2) /= N ) THEN
229

230
          DEALLOCATE( BA )
231
          ALLOCATE( BA(3*SubBand+1,N),stat=istat )
232

233
          IF ( istat /= 0 ) THEN
234
            CALL Fatal( 'BandSolver', 'Memory allocation error.' )
235
          END IF
236

237
       END IF
238

239
       BA = 0.0D0
240
       DO i=1,N
241
         DO j=Rows(i),Rows(i+1)-1
242
           k = i - Cols(j) + 2*Subband + 1
243
           BA(k,Cols(j)) = Values(j)
244
         END DO
245
       END DO
246

247
       CALL SolveBandLapack( N,1,BA,x,Subband,3*Subband+1 )
248

249
     ELSE IF ( A % Format == MATRIX_CRS ) THEN
250

251
       Subband = 0
252
       DO i=1,N
253
         DO j=Rows(i),Diag(i)
254
           Subband = MAX(Subband,ABS(Cols(j)-i))
255
         END DO
256
       END DO
257

258
       IF ( .NOT.ALLOCATED( BA ) ) THEN
259

260
         ALLOCATE( BA(SubBand+1,N),stat=istat )
261

262
         IF ( istat /= 0 ) THEN
263
           CALL Fatal( 'BandSolver', 'Memory allocation error.' )
264
         END IF
265

266
       ELSE IF ( SIZE(BA,1) /= Subband+1 .OR. SIZE(BA,2) /= N ) THEN
267

268
         DEALLOCATE( BA )
269
         ALLOCATE( BA(SubBand+1,N),stat=istat )
270

271
         IF ( istat /= 0 ) THEN
272
           CALL Fatal( 'BandSolver', 'Memory allocation error.' )
273
         END IF
274

275
       END IF
276

277
       BA = 0.0D0
278
       DO i=1,N
279
         DO j=Rows(i),Diag(i)
280
           k = i - Cols(j) + 1
281
           BA(k,Cols(j)) = Values(j)
282
         END DO
283
       END DO
284

285
       CALL SolveSBandLapack( N,1,BA,x,Subband,Subband+1 )
286

287
     ELSE IF ( A % Format == MATRIX_BAND ) THEN
288
       CALL SolveBandLapack( N,1,Values,x,Subband,3*Subband+1 )
289
     ELSE IF ( A % Format == MATRIX_SBAND ) THEN
290
       CALL SolveSBandLapack( N,1,Values,x,Subband,Subband+1 )
291
     END IF
292

293
!------------------------------------------------------------------------------
294
  END SUBROUTINE BandSolver 
295
!------------------------------------------------------------------------------
296

297

298
!------------------------------------------------------------------------------
299
!> Solves a linear system using Umfpack multifrontal direct solver courtesy
300
!> of University of Florida.
301
!------------------------------------------------------------------------------
302
  SUBROUTINE UMFPack_SolveSystem( Solver,A,x,b,Free_Fact )
303
!------------------------------------------------------------------------------
304
    LOGICAL, OPTIONAL :: Free_Fact
305
    TYPE(Matrix_t) :: A
306
    TYPE(Solver_t) :: Solver
307
    REAL(KIND=dp), TARGET :: x(*), b(*)
308

309
    REAL(KIND=dp), POINTER CONTIG :: Values(:)
310
    INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:)
311

312
#include "../config.h"
313
#ifdef HAVE_UMFPACK
314

315
#ifdef ARCH_32_BITS
316
#define CAddrInt c_int32_t
317
#else
318
#define CAddrInt c_int64_t
319
#endif
320
  ! Standard int version
321
  INTERFACE
322
    SUBROUTINE umf4def( control ) &
323
       BIND(C,name='umf4def')
324
       USE, INTRINSIC :: ISO_C_BINDING
325
       REAL(C_DOUBLE) :: control(*)
326
    END SUBROUTINE umf4def
327

328
    SUBROUTINE umf4sym( m,n,rows,cols,values,symbolic,control,iinfo ) &
329
       BIND(C,name='umf4sym')
330
       USE, INTRINSIC :: ISO_C_BINDING
331
       INTEGER(C_INT) :: m,n,rows(*),cols(*)
332
       INTEGER(CAddrInt) ::  symbolic
333
       REAL(C_DOUBLE) :: Values(*), control(*),iinfo(*)
334
    END SUBROUTINE umf4sym
335

336
    SUBROUTINE umf4num( rows,cols,values,symbolic,numeric, control,iinfo ) &
337
       BIND(C,name='umf4num')
338
       USE, INTRINSIC :: ISO_C_BINDING
339
       INTEGER(C_INT) :: rows(*),cols(*)
340
       INTEGER(CAddrInt) ::  numeric, symbolic
341
       REAL(C_DOUBLE) :: Values(*), control(*),iinfo(*)
342
    END SUBROUTINE umf4num
343

344
    SUBROUTINE umf4sol( sys, x, b, numeric, control, iinfo ) &
345
       BIND(C,name='umf4sol')
346
       USE, INTRINSIC :: ISO_C_BINDING
347
       INTEGER(C_INT) :: sys
348
       INTEGER(CAddrInt) :: numeric
349
       REAL(C_DOUBLE) :: x(*), b(*), control(*), iinfo(*)
350
    END SUBROUTINE umf4sol
351

352
    SUBROUTINE umf4fsym(symbolic) &
353
        BIND(C,name='umf4fsym')
354
        USE, INTRINSIC :: ISO_C_BINDING
355
        INTEGER(CAddrInt) :: symbolic
356
    END SUBROUTINE umf4fsym
357

358
    SUBROUTINE umf4fnum(numeric) &
359
        BIND(C,name='umf4fnum')
360
        USE, INTRINSIC :: ISO_C_BINDING
361
        INTEGER(CAddrInt) :: numeric
362
    END SUBROUTINE umf4fnum
363

364
  END INTERFACE
365

366
  ! Long int version
367
  INTERFACE
368
    SUBROUTINE umf4_l_def( control ) &
369
       BIND(C,name='umf4_l_def')
370
       USE, INTRINSIC :: ISO_C_BINDING
371
       REAL(C_DOUBLE) :: control(*)
372
    END SUBROUTINE umf4_l_def
373

374
    SUBROUTINE umf4_l_sym( m,n,rows,cols,values,symbolic,control,iinfo ) &
375
       BIND(C,name='umf4_l_sym')
376
       USE, INTRINSIC :: ISO_C_BINDING
377
       !INTEGER(CAddrInt) ::m,n,rows(*),cols(*) 
378
       INTEGER(UMFPACK_LONG_FORTRAN_TYPE) :: m,n,rows(*),cols(*) !TODO: m,n of are called with AddrInt kind
379
       INTEGER(CAddrInt) ::  symbolic
380
       REAL(C_DOUBLE) :: Values(*), control(*),iinfo(*)
381
    END SUBROUTINE umf4_l_sym
382

383
    SUBROUTINE umf4_l_num( rows,cols,values,symbolic,numeric, control,iinfo ) &
384
       BIND(C,name='umf4_l_num')
385
       USE, INTRINSIC :: ISO_C_BINDING
386
       !INTEGER(CAddrInt) :: rows(*),cols(*)
387
       INTEGER(UMFPACK_LONG_FORTRAN_TYPE) :: rows(*),cols(*)
388
       INTEGER(CAddrInt) ::  numeric, symbolic
389
       REAL(C_DOUBLE) :: Values(*), control(*),iinfo(*)
390
    END SUBROUTINE umf4_l_num
391

392
    SUBROUTINE umf4_l_sol( sys, x, b, numeric, control, iinfo ) &
393
       BIND(C,name='umf4_l_sol')
394
       USE, INTRINSIC :: ISO_C_BINDING
395
       !INTEGER(CAddrInt) :: sys
396
       INTEGER(UMFPACK_LONG_FORTRAN_TYPE) :: sys
397
       INTEGER(CAddrInt) :: numeric
398
       REAL(C_DOUBLE) :: x(*), b(*), control(*), iinfo(*)
399
    END SUBROUTINE umf4_l_sol
400

401
    SUBROUTINE umf4_l_fnum(numeric) &
402
        BIND(C,name='umf4_l_fnum')
403
        USE, INTRINSIC :: ISO_C_BINDING
404
        INTEGER(CAddrInt) :: numeric
405
    END SUBROUTINE umf4_l_fnum
406

407
    SUBROUTINE umf4_l_fsym(symbolic) &
408
        BIND(C,name='umf4_l_fsym')
409
        USE, INTRINSIC :: ISO_C_BINDING
410
        INTEGER(CAddrInt) :: symbolic
411
    END SUBROUTINE umf4_l_fsym
412
  END INTERFACE
413

414
  INTEGER :: i, n, status, sys
415
  REAL(KIND=dp) :: iInfo(90), Control(20)
416
  INTEGER(KIND=AddrInt) :: symbolic, zero=0
417
  INTEGER(KIND=UMFPACK_LONG_FORTRAN_TYPE) :: ln, lsys
418
  INTEGER(KIND=UMFPACK_LONG_FORTRAN_TYPE), ALLOCATABLE :: LRows(:), LCols(:)
419

420
  SAVE iInfo, Control
421
 
422
  LOGICAL :: Factorize, FreeFactorize, stat, BigMode
423

424
  IF ( PRESENT(Free_Fact) ) THEN
425
    IF ( Free_Fact ) THEN
426
      IF ( A % UMFPack_Numeric/=0 ) THEN
427
        CALL umf4fnum(A % UMFPack_Numeric)
428
        A % UMFPack_Numeric = 0
429
      END IF
430
      RETURN
431
    END IF
432
  END IF
433

434
  BigMode = ListGetString( Solver % Values, &
435
      'Linear System Direct Method' ) == 'big umfpack'
436

437
  Factorize = ListGetLogical( Solver % Values, &
438
     'Linear System Refactorize', stat )
439
  IF ( .NOT. stat ) Factorize = .TRUE.
440

441
  n = A % NumberofRows
442
  Rows => A % Rows
443
  Cols => A % Cols
444
  Diag => A % Diag
445
  Values => A % Values
446

447
  IF ( Factorize .OR. A% UmfPack_Numeric==0 ) THEN
448
    IF ( A % UMFPack_Numeric /= 0 ) THEN
449
      IF( BigMode ) THEN
450
        CALL umf4_l_fnum( A % UMFPack_Numeric )
451
      ELSE
452
        CALL umf4fnum( A % UMFPack_Numeric )
453
      END IF
454
      A % UMFPack_Numeric = 0
455
    END IF
456

457
    IF ( BigMode ) THEN
458
      ALLOCATE( LRows(SIZE(Rows)), LCols(SIZE(Cols)) )
459

460
      DO i=1,n+1
461
        LRows(i) = Rows(i)-1
462
      END DO
463

464
      DO i=1,SIZE(Cols)
465
        LCols(i) = Cols(i)-1
466
      END DO
467

468
      ln = n ! TODO: Kludge: ln is AddrInt and n is regular INTEGER
469
      CALL umf4_l_def( Control )
470
      CALL umf4_l_sym( ln,ln, LRows, LCols, Values, Symbolic, Control, iInfo )
471
    ELSE
472
      Rows = Rows-1
473
      Cols = Cols-1
474
      CALL umf4def( Control )
475
      CALL umf4sym( n,n, Rows, Cols, Values, Symbolic, Control, iInfo )
476
    END IF
477

478
    IF (iInfo(1)<0) THEN
479
      PRINT *, 'Error occurred in umf4sym: ', iInfo(1)
480
      STOP EXIT_ERROR
481
    END IF
482

483
    IF ( BigMode ) THEN
484
      CALL umf4_l_num(LRows, LCols, Values, Symbolic, A % UMFPack_Numeric, Control, iInfo )
485
    ELSE
486
      CALL umf4num( Rows, Cols, Values, Symbolic, A % UMFPack_Numeric, Control, iInfo )
487
    END IF
488

489
    IF (iinfo(1)<0) THEN
490
      PRINT*, 'Error occurred in umf4num: ', iinfo(1)
491
      STOP EXIT_ERROR
492
    ENDIF
493

494
    IF ( BigMode ) THEN
495
      DEALLOCATE( LRows, LCols )
496
      CALL umf4_l_fsym( Symbolic )
497
    ELSE
498
      A % Rows = A % Rows+1
499
      A % Cols = A % Cols+1
500
      CALL umf4fsym( Symbolic )
501
    END IF
502
  END IF
503

504
  IF ( BigMode ) THEN
505
    lsys = 2
506
    CALL umf4_l_sol( lsys, x, b, A % UMFPack_Numeric, Control, iInfo )
507
  ELSE
508
    sys = 2
509
    CALL umf4sol( sys, x, b, A % UMFPack_Numeric, Control, iInfo )
510
  END IF
511

512
  IF (iinfo(1)<0) THEN
513
    PRINT*, 'Error occurred in umf4sol: ', iinfo(1)
514
    STOP EXIT_ERROR
515
  END IF
516
 
517
  FreeFactorize = ListGetLogical( Solver % Values, &
518
      'Linear System Free Factorization', stat )
519
  IF ( .NOT. stat ) FreeFactorize = .TRUE.
520

521
  IF ( Factorize .AND. FreeFactorize ) THEN
522
    IF ( BigMode ) THEN
523
      CALL umf4_l_fnum(A % UMFPack_Numeric)
524
    ELSE
525
      CALL umf4fnum(A % UMFPack_Numeric)
526
    END IF
527
    A % UMFPack_Numeric = 0
528
  END IF
529
#else
530
   CALL Fatal( 'UMFPack_SolveSystem', 'UMFPACK Solver has not been installed.' )
531
#endif
532
!------------------------------------------------------------------------------
533
  END SUBROUTINE UMFPack_SolveSystem
534
!------------------------------------------------------------------------------
535

536

537
!------------------------------------------------------------------------------
538
!> Solves a linear system using Cholmod multifrontal direct solver courtesy
539
!> of University of Florida.
540
!------------------------------------------------------------------------------
541
  SUBROUTINE Cholmod_SolveSystem( Solver,A,x,b,Free_fact)
542
!------------------------------------------------------------------------------
543
  LOGICAL, OPTIONAL :: Free_Fact
544
  TYPE(Matrix_t) :: A
545
  TYPE(Solver_t) :: Solver
546
  REAL(KIND=dp) :: x(*), b(*)
547

548
  INTERFACE
549
     SUBROUTINE cholmod_ffree(chol) BIND(c,NAME="cholmod_ffree")
550
       USE Types
551
       INTEGER(KIND=AddrInt) :: chol
552
     END SUBROUTINE cholmod_ffree
553

554
     FUNCTION cholmod_ffactorize(n,rows,cols,vals,cmplx) RESULT(chol) BIND(c,NAME="cholmod_ffactorize")
555
        USE Types
556
        INTEGER :: n, cmplx, Rows(*), Cols(*)
557
        REAL(KIND=dp) :: Vals(*)
558
        INTEGER(KIND=dp) :: chol
559
     END FUNCTION cholmod_ffactorize
560

561
     SUBROUTINE cholmod_fsolve(chol, n, x,b) BIND(c,NAME="cholmod_fsolve")
562
        USE Types
563
        REAL(KIND=dp) :: x(*), b(*)
564
        INTEGER :: n
565
        INTEGER(KIND=dp) :: chol
566
     END SUBROUTINE cholmod_fsolve
567
  END INTERFACE
568

569
  LOGICAL :: Factorize, FreeFactorize, Found
570
  INTEGER :: i
571
  REAL(KIND=dp), POINTER CONTIG :: Vals(:)
572
  INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:)
573

574
#ifdef HAVE_CHOLMOD
575
  IF ( PRESENT(Free_Fact) ) THEN
576
    IF ( Free_Fact ) THEN
577
      IF ( A % Cholmod/=0 ) THEN
578
        CALL cholmod_ffree(A % cholmod)
579
        A % cholmod = 0
580
      END IF
581
      RETURN
582
    END IF
583
  END IF
584

585
  Factorize = ListGetLogical( Solver % Values, &
586
     'Linear System Refactorize', Found )
587
  IF ( .NOT. Found ) Factorize = .TRUE.
588

589
  IF ( Factorize .OR. A% cholmod==0 ) THEN
590
    IF ( A % cholmod/=0 ) THEN
591
      CALL cholmod_ffree(A % cholmod)
592
      A % cholmod = 0
593
    END IF
594

595
    Rows => A % Rows
596
    Cols => A % Cols
597
    Vals => A % Values
598

599
    Rows=Rows-1; Cols=Cols-1 ! c numbering
600
    i=0
601
    IF(A % Complex) i=1
602
    A % Cholmod=cholmod_ffactorize(A % NumberOfRows, Rows, Cols, Vals, i)
603
    Rows=Rows+1; Cols=Cols+1 ! fortran numbering
604
  END IF
605

606
  CALL cholmod_fsolve(A % cholmod, A % NumberOfRows, x, b);
607

608
  FreeFactorize = ListGetLogical( Solver % Values, &
609
      'Linear System Free Factorization', Found )
610
  IF ( .NOT. Found ) FreeFactorize = .TRUE.
611

612
  IF ( Factorize .AND. FreeFactorize ) THEN
613
    CALL cholmod_ffree(A % cholmod)
614
    A % cholmod = 0
615
  END IF
616
#else
617
   CALL Fatal( 'Cholmod_SolveSystem', 'Cholmod Solver has not been installed.' )
618
#endif
619
!------------------------------------------------------------------------------
620
  END SUBROUTINE Cholmod_SolveSystem
621
!------------------------------------------------------------------------------
622

623

624
!------------------------------------------------------------------------------
625
!> Solves a linear system using SuiteSparseQR multifrontal direct solver courtesy
626
!> of University of Florida.
627
!------------------------------------------------------------------------------
628
  SUBROUTINE SPQR_SolveSystem( Solver,A,x,b,Free_fact)
629
!------------------------------------------------------------------------------
630
  LOGICAL, OPTIONAL :: Free_Fact
631
  TYPE(Matrix_t) :: A
632
  TYPE(Solver_t) :: Solver
633
  REAL(KIND=dp) :: x(*), b(*)
634

635
  INTEGER :: i
636
  LOGICAL :: Factorize, FreeFactorize, Found
637

638
  REAL(KIND=dp), POINTER CONTIG :: Vals(:)
639
  INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:)
640

641
  INTERFACE
642
     FUNCTION spqr_ffree(chol) RESULT(stat) BIND(c,NAME="spqr_ffree")
643
       USE Types
644
       INTEGER :: stat
645
       INTEGER(KIND=AddrInt) :: chol
646
     END FUNCTION spqr_ffree
647

648
     FUNCTION spqr_ffactorize(n,rows,cols,vals) RESULT(chol) BIND(c,NAME="spqr_ffactorize")
649
       USE Types
650
       INTEGER :: n, rows(*), cols(*)
651
       REAL(KIND=dp) :: vals(*)
652
       INTEGER(KIND=AddrInt) :: chol
653
     END FUNCTION spqr_ffactorize
654

655
     SUBROUTINE spqr_fsolve(chol, n, x,b) BIND(c,NAME="spqr_fsolve")
656
       USE Types
657
       REAL(KIND=dp) :: x(*), b(*)
658
       INTEGER :: n
659
       INTEGER(KIND=AddrInt) :: chol
660
     END SUBROUTINE spqr_fsolve
661
  END INTERFACE
662

663
#ifdef HAVE_CHOLMOD
664
  IF ( PRESENT(Free_Fact) ) THEN
665
    IF ( Free_Fact ) THEN
666
      IF ( A % Cholmod/=0 ) THEN
667
        IF(spqr_ffree(A % cholmod)==0) A % Cholmod=0
668
      END IF
669
      RETURN
670
    END IF
671
  END IF
672

673
  Factorize = ListGetLogical( Solver % Values, &
674
     'Linear System Refactorize', Found )
675
  IF ( .NOT. Found ) Factorize = .TRUE.
676

677
  IF ( Factorize .OR. A% cholmod==0 ) THEN
678
    IF ( A % cholmod/=0 ) THEN
679
      i=spqr_ffree(A % cholmod)
680
      A % cholmod = 0
681
    END IF
682

683
    Rows => A % Rows
684
    Cols => A % Cols
685
    Vals => A % Values
686

687
    Rows=Rows-1; Cols=Cols-1 ! c numbering
688
    A % Cholmod=spqr_ffactorize(A % NumberOfRows, Rows, Cols, Vals)
689
    Rows=Rows+1; Cols=Cols+1 ! fortran numbering
690
  END IF
691

692
  CALL spqr_fsolve(A % cholmod, A % NumberOfRows, x, b);
693

694
  FreeFactorize = ListGetLogical( Solver % Values, &
695
      'Linear System Free Factorization', Found )
696
  IF ( .NOT. Found ) FreeFactorize = .TRUE.
697

698
  IF ( Factorize .AND. FreeFactorize ) THEN
699
    i=spqr_ffree(A % cholmod)
700
    A % cholmod = 0
701
  END IF
702
#else
703
   CALL Fatal( 'SPQR_SolveSystem', 'SPQR Solver has not been installed.' )
704
#endif
705
!------------------------------------------------------------------------------
706
  END SUBROUTINE SPQR_SolveSystem
707
!------------------------------------------------------------------------------
708

709
!------------------------------------------------------------------------------
710
 SUBROUTINE FreeMumpsFactorizations(A)
711
!------------------------------------------------------------------------------
712
    TYPE(Matrix_t) :: A
713

714
!------------------------------------------------------------------------------
715
#ifdef HAVE_MUMPS
716
    IF(ASSOCIATED(A % SMumpsID)) THEN
717
      DEALLOCATE( A % SMumpsID % irn_loc, &
718
         A % SMumpsID % jcn_loc, A % SMumpsID % Rhs,  &
719
           A % SMumpsID % isol_loc, A % SMumpsID % sol_loc)
720

721
      DEALLOCATE( A % SMumpsID % A_loc )
722
      IF (ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
723
      A % Gorder=>Null()
724

725
      A % SMumpsID % job = -2
726
      CALL SMumps(A % SMumpsID)
727
      DEALLOCATE(A % SMumpsID)
728
      A % SMumpsId => NULL()
729
    END IF
730

731
    IF(ASSOCIATED(A % CMumpsID)) THEN
732
      DEALLOCATE( A % CMumpsID % irn_loc, &
733
         A % CMumpsID % jcn_loc, A % CMumpsID % Rhs,  &
734
           A % CMumpsID % isol_loc, A % CMumpsID % sol_loc)
735

736
      DEALLOCATE( A % CMumpsID % A_loc )
737
      IF (ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
738
      A % Gorder=>Null()
739

740
      A % CMumpsID % job = -2
741
      CALL CMumps(A % CMumpsID)
742
      DEALLOCATE(A % CMumpsID)
743
      A % CMumpsId => NULL()
744
    END IF
745

746
    IF(ASSOCIATED(A % MumpsID)) THEN
747
      DEALLOCATE( A % MumpsID % irn_loc, &
748
         A % MumpsID % jcn_loc, A % MumpsID % Rhs,  &
749
           A % MumpsID % isol_loc, A % MumpsID % sol_loc)
750

751
      IF(.NOT.ASSOCIATED(A % MumpsID % A_loc, A % Values)) DEALLOCATE( A % MumpsID % A_loc )
752
      IF (ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
753
      A % Gorder=>Null()
754

755
      A % MumpsID % job = -2
756
      CALL DMumps(A % MumpsID)
757
      DEALLOCATE(A % MumpsID)
758
      A % MumpsId => NULL()
759
    END IF
760

761
    IF(ASSOCIATED(A % ZMumpsID)) THEN
762
      DEALLOCATE( A % ZMumpsID % irn_loc, &
763
         A % ZMumpsID % jcn_loc, A % ZMumpsID % Rhs,  &
764
           A % ZMumpsID % isol_loc, A % ZMumpsID % sol_loc)
765

766
      DEALLOCATE( A % ZMumpsID % A_loc )
767

768
      IF (ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
769
      A % Gorder=>Null()
770

771
      A % ZMumpsID % job = -2
772
      CALL ZMumps(A % ZMumpsID)
773
      DEALLOCATE(A % ZMumpsID)
774
      A % ZMumpsId => NULL()
775
    END IF
776
#endif
777
!------------------------------------------------------------------------------
778
 END SUBROUTINE FreeMumpsFactorizations
779
!------------------------------------------------------------------------------
780

781

782
!------------------------------------------------------------------------------
783
!> Solves a linear system using MUMPS direct solver. This is a legacy solver
784
!> with complicated dependencies. Single precision version.
785
!------------------------------------------------------------------------------
786
  SUBROUTINE SMumps_SolveSystem( Solver,A,x,b )
787
!------------------------------------------------------------------------------
788
#ifdef HAVE_MUMPS
789
#  if defined(ELMER_HAVE_MPI_MODULE)
790
  USE mpi
791
#  endif
792
#endif
793

794
  TYPE(Matrix_t) :: A
795
  TYPE(Solver_t) :: Solver
796
  REAL(KIND=dp), TARGET :: x(*), b(*)
797

798
#ifdef HAVE_MUMPS
799
#  if defined(ELMER_HAVE_MPIF_HEADER)
800
  INCLUDE 'mpif.h'
801
#  endif
802

803
  INTEGER, ALLOCATABLE :: Owner(:)
804
  INTEGER :: i,j,n,ip,ierr,icntlft,nzloc
805
  LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled
806

807
  INTEGER, ALLOCATABLE :: memb(:)
808
  INTEGER :: Comm_active, Group_active, Group_world
809

810
  REAL, ALLOCATABLE :: dbuf(:)
811

812
  Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat )
813
  IF ( .NOT. stat ) Factorize = .TRUE.
814

815
  IF ( Factorize .OR. .NOT.ASSOCIATED(A % SMumpsID) ) THEN
816
    CALL FreeMumpsFactorizations(A)
817
    ALLOCATE(A % SMumpsID)
818

819
    A % SMumpsID % Comm = A % Comm
820
    A % SMumpsID % par  =  1
821
    A % SMumpsID % job  = -1
822
    A % SMumpsID % Keep =  0
823

824
    matsym = ListGetLogical( Solver % Values, 'Linear System Symmetric', stat)
825
    matspd = ListGetLogical( Solver % Values, 'Linear System Positive Definite', stat)
826

827
    ! force unsymmetric mode when "row equilibration" is used
828
    scaled = ListGetLogical( Solver % Values, 'Linear System Scaling', stat)
829
    IF(.NOT.stat) scaled = .TRUE.
830
    IF(scaled) THEN
831
      IF(ListGetLogical( Solver % Values, 'Linear System Row Equilibration',stat)) matsym=.FALSE.
832
    END IF
833

834
    IF(matsym) THEN
835
      IF ( matspd) THEN
836
        A % MumpsID % sym = 1
837
      ELSE
838
        A % SMumpsID % sym = 0 ! 2=symmetric, but unsymmetric solver seems faster, at least in a few
839
                              ! simple cases...  more testing needed...
840
      END IF
841
    ELSE
842
      A % SMumpsID % sym = 0
843
    END IF
844

845
    CALL SMumps(A % SMumpsID)
846

847
    IF(ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
848

849
    IF(ASSOCIATED(A % ParallelInfo)) THEN
850
      n = SIZE(A % ParallelInfo % GlobalDOFs)
851

852
      ALLOCATE( A % Gorder(n), Owner(n) )
853
      CALL ContinuousNumbering( A % ParallelInfo, A % Perm, A % Gorder, Owner )
854

855
      CALL MPI_ALLREDUCE( SUM(Owner), A % SMumpsID % n, &
856
         1, MPI_INTEGER, MPI_SUM, A % SMumpsID % Comm, ierr )
857
      DEALLOCATE(Owner)
858
    ELSE
859
      CALL MPI_ALLREDUCE( A % NumberOfRows, A % SMumpsId % n, &
860
            1, MPI_INTEGER, MPI_MAX, A % Comm, ierr )
861

862
      ALLOCATE(A % Gorder(A % NumberOFrows))
863
      DO i=1,A % NumberOFRows
864
        A % Gorder(i) = i
865
      END DO
866
    END IF
867

868
   ! Set matrix for Mumps (unsymmetric case)
869
    IF (A % SmumpsID % sym == 0) THEN
870
      A % SMumpsID % nz_loc = A % Rows(A % NumberOfRows+1)-1
871

872
      ALLOCATE( A % SMumpsID % irn_loc(A % SMumpsID % nz_loc) )
873
      ALLOCATE( A % SMumpsID % a_loc(A % SMumpsId % nz_loc) )
874
      ALLOCATE( A % SMumpsID % jcn_loc(A % SMumpsId % nz_loc) )
875

876
      nzloc = 0
877
      DO i=1,A % NumberOfRows
878
        ip = A % Gorder(i)
879
        DO j=A % Rows(i),A % Rows(i+1)-1
880
          nzloc = nzloc + 1
881
          A % SMumpsID % irn_loc(nzloc) = ip
882
          A % SMumpsID % a_loc(nzloc)   = A % Values(j)
883
          A % SMumpsID % jcn_loc(nzloc) = A % Gorder(A % Cols(j))
884
        END DO
885
      END DO
886
    ELSE
887
      ! Set matrix for Mumps (symmetric case)
888
      nzloc = 0
889
      DO i=1,A % NumberOfRows
890
        ! Only output lower triangular part to Mumps
891
        DO j=A % Rows(i),A % Diag(i)
892
          nzloc = nzloc + 1
893
        END DO
894
      END DO
895

896
      A % SMumpsID % nz_loc = nzloc
897

898
      ALLOCATE( A % SMumpsID % irn_loc(A % SMumpsID % nz_loc) )
899
      ALLOCATE( A % SMumpsID % jcn_loc(A % SMumpsId % nz_loc) )
900
      ALLOCATE( A % SMumpsID % A_loc(A % SMumpsId % nz_loc) )
901

902
      nzloc = 0
903
      DO i=1,A % NumberOfRows
904
        ! Only output lower triangular part to Mumps
905
        ip = A % Gorder(i)
906
        DO j=A % Rows(i),A % Diag(i)
907
          nzloc = nzloc + 1
908
          A % SmumpsID % IRN_loc(nzloc) = ip
909
          A % SmumpsID % A_loc(nzloc) = A % Values(j)
910
          A % SmumpsID % JCN_loc(nzloc) = A % Gorder(A % Cols(j))
911
        END DO
912
      END DO
913
    END IF
914

915

916
    ALLOCATE(A % SMumpsID % rhs(A % SMumpsId % n))
917

918
    ! Tune verbosity of MUMPS.
919
    i = 0
920
    IF(InfoActive(20)) i = 1
921
    A % SMumpsID % icntl(2) = i ! suppress printing of diagnostics and warnings
922
    A % SMumpsID % icntl(3) = i ! suppress statistics
923

924
    A % SMumpsID % icntl(4) = 1 ! the same as the two above, but doesn't seem to work.
925
    A % SMumpsID % icntl(5) = 0 ! matrix format 'assembled'
926

927
    icntlft = ListGetInteger(Solver % Values, &
928
          'mumps percentage increase working space', stat)
929
    IF (stat) THEN
930
       A % SMumpsID % icntl(14) = icntlft
931
    END IF
932
    A % SMumpsID % icntl(18) = 3 ! 'distributed' matrix 
933
    A % SMumpsID % icntl(21) = 1 ! 'distributed' solution phase
934

935
    A % SMumpsID % job = 4
936
    CALL SMumps(A % SMumpsID)
937
    CALL Flush(6)
938

939
    A % SMumpsID % lsol_loc = A % Smumpsid % info(23)
940
    ALLOCATE(A % SMumpsID % sol_loc(A % SMumpsId % lsol_loc))
941
    ALLOCATE(A % SMumpsID % isol_loc(A % SMumpsId % lsol_loc))
942
  END IF
943

944
 ! sum the rhs from all procs. Could be done for neighbours only (i guess):
945
 ! ------------------------------------------------------------------------
946
  A % SMumpsID % RHS = 0
947
  DO i=1,A % NumberOfRows
948
    ip = A % Gorder(i)
949
    A % SMumpsId % RHS(ip) = b(i)
950
  END DO
951

952
  ALLOCATE( dbuf(A % SMumpsID % n) )
953
  dbuf = A % SMumpsId % RHS
954
  CALL MPI_ALLREDUCE( dbuf, A % SMumpsID % RHS, &
955
    A % SMumpsID % n, MPI_REAL, MPI_SUM, A % SMumpsID % Comm, ierr )
956

957
 ! Solution:
958
 ! ---------
959
  A % SMumpsID % job = 3
960
  CALL SMumps(A % SMumpsID)
961

962
 ! Distribute the solution to all:
963
 ! -------------------------------
964
  A % SMumpsId % Rhs = 0
965
  DO i=1,A % SMumpsID % lsol_loc
966
    A % SMumpsID % RHS(A % SMumpsID % isol_loc(i)) = A % SMumpsID % sol_loc(i)
967
  END DO
968
  dbuf = A % SMumpsId % RHS
969
  CALL MPI_ALLREDUCE( dbuf, A % SMumpsID % RHS, &
970
    A % SMumpsID % n, MPI_REAL, MPI_SUM,A %  SMumpsID % Comm, ierr )
971

972
  DEALLOCATE(dbuf)
973

974
 ! Select the values which belong to us:
975
 ! -------------------------------------
976
  DO i=1,A % NumberOfRows
977
    ip = A % Gorder(i)
978
    x(i) = A % SMumpsId % RHS(ip)
979
  END DO
980

981
  FreeFactorize = ListGetLogical( Solver % Values, 'Linear System Free Factorization', stat )
982
  IF ( .NOT. stat ) FreeFactorize = .TRUE.
983

984
  IF ( Factorize .AND. FreeFactorize ) CALL FreeMumpsFactorizations(A)
985
#else
986
   CALL Fatal( 'Mumps_SolveSystem', 'MUMPS Solver has not been installed.' )
987
#endif
988
!------------------------------------------------------------------------------
989
  END SUBROUTINE SMumps_SolveSystem
990
!------------------------------------------------------------------------------
991

992

993
!------------------------------------------------------------------------------
994
!> Solves a linear system using MUMPS direct solver. This is a legacy solver
995
!> with complicated dependencies. Single precision complex version.
996
!------------------------------------------------------------------------------
997
  SUBROUTINE CMumps_SolveSystem( Solver,A,x,b )
998
!------------------------------------------------------------------------------
999
#ifdef HAVE_MUMPS
1000
#  if defined(ELMER_HAVE_MPI_MODULE)
1001
  USE mpi
1002
#  endif
1003
#endif
1004

1005
  TYPE(Matrix_t) :: A
1006
  TYPE(Solver_t) :: Solver
1007
  REAL(KIND=dp), TARGET :: x(*), b(*)
1008

1009
#ifdef HAVE_MUMPS
1010
#  if defined(ELMER_HAVE_MPIF_HEADER)
1011
  INCLUDE 'mpif.h'
1012
#  endif
1013

1014
  INTEGER, ALLOCATABLE :: Owner(:)
1015
  INTEGER :: i,j,n,ip,ierr,icntlft,nzloc
1016
  LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled
1017

1018
  INTEGER, ALLOCATABLE :: memb(:)
1019
  INTEGER :: Comm_active, Group_active, Group_world
1020

1021
  COMPLEX, ALLOCATABLE :: dbuf(:)
1022

1023
  Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat )
1024
  IF ( .NOT. stat ) Factorize = .TRUE.
1025

1026
  IF ( Factorize .OR. .NOT.ASSOCIATED(A % CMumpsID) ) THEN
1027
    CALL FreeMumpsFactorizations(A)
1028
    ALLOCATE(A % CMumpsID)
1029

1030
    A % CMumpsID % Comm = A % Comm
1031
    A % CMumpsID % par  =  1
1032
    A % CMumpsID % job  = -1
1033
    A % CMumpsID % Keep =  0
1034

1035
!   matsym = ListGetLogical( Solver % Values, 'Linear System Symmetric', stat)
1036
!   matspd = ListGetLogical( Solver % Values, 'Linear System Positive Definite', stat)
1037
    matsym = .FALSE.
1038
    matspd = .FALSE.
1039

1040
    ! force unsymmetric mode when "row equilibration" is used
1041
    scaled = ListGetLogical( Solver % Values, 'Linear System Scaling', stat)
1042
    IF(.NOT.stat) scaled = .TRUE.
1043
    IF(scaled) THEN
1044
      IF(ListGetLogical( Solver % Values, 'Linear System Row Equilibration',stat)) matsym=.FALSE.
1045
    END IF
1046

1047
!   IF(matsym) THEN
1048
!     IF ( matspd) THEN
1049
!       A % CMumpsID % sym = 1
1050
!     ELSE
1051
!       A % CMumpsID % sym = 0 ! 2=symmetric, but unsymmetric solver seems faster, at least in a few
1052
!                             ! simple cases...  more testing needed...
1053
!     END IF
1054
!   ELSE
1055
!     A % CMumpsID % sym = 0
1056
!   END IF
1057
     A % CMumpsID % sym = 0
1058

1059
    CALL CMumps(A % CMumpsID)
1060

1061
    IF(ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
1062

1063
    IF(ASSOCIATED(A % ParallelInfo)) THEN
1064
      n = SIZE(A % ParallelInfo % GlobalDOFs)
1065
      ALLOCATE( A % Gorder(n), Owner(n) )
1066
      CALL ContinuousNumbering( A % ParallelInfo, A % Perm, A % Gorder, Owner )
1067

1068
      CALL MPI_ALLREDUCE( SUM(Owner)/2, A % CMumpsID % n, &
1069
        1, MPI_INTEGER, MPI_SUM, A % CMumpsID % Comm, ierr )
1070
      DEALLOCATE(Owner)
1071
    ELSE
1072
      CALL MPI_ALLREDUCE( A % NumberOfRows/2, A % CMumpsId % n, &
1073
            1, MPI_INTEGER, MPI_MAX, A % Comm, ierr )
1074

1075
      ALLOCATE(A % Gorder(A % NumberOFrows))
1076
      DO i=1,A % NumberOFRows
1077
        A % Gorder(i) = i
1078
      END DO
1079
    END IF
1080

1081
   ! Set matrix for Mumps (unsymmetric case)
1082
    IF (A % CmumpsID % sym == 0) THEN
1083
      A % CMumpsID % nz_loc = (A % Rows(A % NumberOfRows+1)-1)/4
1084

1085
      ALLOCATE( A % CMumpsID % irn_loc(A % CMumpsID % nz_loc) )
1086
      ALLOCATE( A % CMumpsID % a_loc(A % CMumpsId % nz_loc) )
1087
      ALLOCATE( A % CMumpsID % jcn_loc(A % CMumpsId % nz_loc) )
1088

1089
      nzloc = 0
1090
      DO i=1,A % NumberOfRows,2
1091
        ip = (A % Gorder(i)-1)/2+1
1092
        DO j=A % Rows(i),A % Rows(i+1)-1,2
1093
          nzloc = nzloc + 1
1094
          A % CMumpsID % irn_loc(nzloc) = ip
1095
          A % CMumpsID % jcn_loc(nzloc) = (A % Gorder(A % Cols(j))-1)/2+1
1096
          A % CMumpsID % a_loc(nzloc)   = CMPLX( A % Values(j), -A % Values(j+1) )
1097
        END DO
1098
      END DO
1099
    ELSE
1100
      ! Set matrix for Mumps (symmetric case)
1101
      nzloc = 0
1102
      DO i=1,A % NumberOfRows
1103
        ! Only output lower triangular part to Mumps
1104
        DO j=A % Rows(i),A % Diag(i)
1105
          nzloc = nzloc + 1
1106
        END DO
1107
      END DO
1108

1109
      A % CMumpsID % nz_loc = nzloc
1110

1111
      ALLOCATE( A % CMumpsID % irn_loc(A % CMumpsID % nz_loc) )
1112
      ALLOCATE( A % CMumpsID % jcn_loc(A % CMumpsId % nz_loc) )
1113
      ALLOCATE( A % CMumpsID % A_loc(A % CMumpsId % nz_loc) )
1114

1115
      nzloc = 0
1116
      DO i=1,A % NumberOfRows,2
1117
        ! Only output lower triangular part to Mumps
1118
        ip = A % Gorder(i)
1119
        DO j=A % Rows(i),A % Diag(i),2
1120
          nzloc = nzloc + 1
1121
          A % CmumpsID % IRN_loc(nzloc) = ip
1122
          A % CmumpsID % A_loc(nzloc) = A % Values(j)
1123
          A % CmumpsID % JCN_loc(nzloc) = A % Gorder(A % Cols(j))
1124
        END DO
1125
      END DO
1126
    END IF
1127

1128

1129
    ALLOCATE(A % CMumpsID % rhs(A % CMumpsId % n))
1130

1131
    ! Tune verbosity of MUMPS.
1132
    i = 0
1133
    IF(InfoActive(20)) i = 1
1134
    A % CMumpsID % icntl(2) = i ! suppress printing of diagnostics and warnings
1135
    A % CMumpsID % icntl(3) = i ! suppress statistics
1136

1137
    A % CMumpsID % icntl(4) = 1 ! the same as the two above, but doesn't seem to work.
1138
    A % CMumpsID % icntl(5) = 0 ! matrix format 'assembled'
1139

1140
    icntlft = ListGetInteger(Solver % Values, &
1141
          'mumps percentage increase working space', stat)
1142
    IF (stat) THEN
1143
       A % CMumpsID % icntl(14) = icntlft
1144
    END IF
1145
    A % CMumpsID % icntl(18) = 3 ! 'distributed' matrix 
1146
    A % CMumpsID % icntl(21) = 1 ! 'distributed' solution phase
1147

1148
    A % CMumpsID % job = 4
1149
    CALL CMumps(A % CMumpsID)
1150
    CALL Flush(6)
1151

1152
    A % CMumpsID % lsol_loc = A % CMumpsid % info(23)
1153
    ALLOCATE(A % CMumpsID % sol_loc(A % CMumpsId % lsol_loc))
1154
    ALLOCATE(A % CMumpsID % isol_loc(A % CMumpsId % lsol_loc))
1155
  END IF
1156

1157
 ! sum the rhs from all procs. Could be done
1158
 ! for neighbours only (i guess):
1159
 ! ------------------------------------------
1160
  A % CMumpsID % RHS = 0
1161
  DO i=1,A % NumberOfRows,2
1162
    ip = (A % Gorder(i)-1)/2+1
1163
    A % CMumpsId % RHS(ip) = CMPLX( b(i), b(i+1) )
1164
  END DO
1165

1166
  ALLOCATE( dbuf(A % CMumpsID % n) )
1167
  dbuf = A % CMumpsId % RHS
1168
  CALL MPI_ALLREDUCE( dbuf, A % CMumpsID % RHS, &
1169
    A % CMumpsID % n, MPI_COMPLEX, MPI_SUM, A % CMumpsID % Comm, ierr )
1170

1171
 ! Solution:
1172
 ! ---------
1173
  A % CMumpsID % job = 3
1174
  CALL CMumps(A % CMumpsID)
1175

1176
 ! Distribute the solution to all:
1177
 ! -------------------------------
1178
  A % CMumpsId % Rhs = 0
1179
  DO i=1,A % CMumpsID % lsol_loc
1180
    A % CMumpsID % RHS(A % CMumpsID % isol_loc(i)) = A % CMumpsID % sol_loc(i)
1181
  END DO
1182
  dbuf = A % CMumpsId % RHS
1183
  CALL MPI_ALLREDUCE( dbuf, A % CMumpsID % RHS, &
1184
    A % CMumpsID % N, MPI_COMPLEX, MPI_SUM, A %  CMumpsID % Comm, ierr )
1185

1186
  DEALLOCATE(dbuf)
1187

1188
 ! Select the values which belong to us:
1189
 ! -------------------------------------
1190
  DO i=1,A % NumberOfRows,2
1191
    ip = (A % Gorder(i)-1)/2+1
1192
    x(i)   = REAL( A % CMumpsId % RHS(ip) )
1193
    x(i+1) = AIMAG( A % CMumpsId % RHS(ip) )
1194
  END DO
1195

1196
  FreeFactorize = ListGetLogical( Solver % Values, 'Linear System Free Factorization', stat )
1197
  IF ( .NOT. stat ) FreeFactorize = .TRUE.
1198
  IF ( Factorize .AND. FreeFactorize ) CALL FreeMumpsFactorizations(A)
1199

1200
#else
1201
   CALL Fatal( 'Mumps_SolveSystem', 'MUMPS Solver has not been installed.' )
1202
#endif
1203
!------------------------------------------------------------------------------
1204
  END SUBROUTINE CMumps_SolveSystem
1205
!------------------------------------------------------------------------------
1206

1207
!------------------------------------------------------------------------------
1208
!> Solves a linear system using MUMPS direct solver. This is a legacy solver
1209
!> with complicated dependencies. This is only available in parallel. 
1210
!------------------------------------------------------------------------------
1211
  SUBROUTINE Mumps_SolveSystem( Solver,A,x,b )
1212
!------------------------------------------------------------------------------
1213
#ifdef HAVE_MUMPS
1214
#  if defined(ELMER_HAVE_MPI_MODULE)
1215
  USE mpi
1216
#  endif
1217
#endif
1218

1219
  TYPE(Matrix_t) :: A
1220
  TYPE(Solver_t) :: Solver
1221
  REAL(KIND=dp), TARGET :: x(*), b(*)
1222

1223
#ifdef HAVE_MUMPS
1224
#  if defined(ELMER_HAVE_MPIF_HEADER)
1225
  INCLUDE 'mpif.h'
1226
#  endif
1227

1228
  INTEGER, ALLOCATABLE :: Owner(:)
1229
  INTEGER :: i,j,n,ip,ierr,icntlft,nzloc
1230
  LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled
1231

1232
  INTEGER, ALLOCATABLE :: memb(:)
1233
  INTEGER :: Comm_active, Group_active, Group_world
1234

1235
  REAL(KIND=dp), ALLOCATABLE :: dbuf(:)
1236

1237

1238
  Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat )
1239
  IF ( .NOT. stat ) Factorize = .TRUE.
1240

1241
  IF ( Factorize .OR. .NOT.ASSOCIATED(A % MumpsID) ) THEN
1242
    CALL FreeMumpsFactorizations(A)
1243
    ALLOCATE(A % MumpsID)
1244

1245
    A % MumpsID % Comm = A % Comm
1246
    A % MumpsID % par  =  1
1247
    A % MumpsID % job  = -1
1248
    A % MumpsID % Keep =  0
1249

1250
    matsym = ListGetLogical( Solver % Values, 'Linear System Symmetric', stat)
1251
    matspd = ListGetLogical( Solver % Values, 'Linear System Positive Definite', stat)
1252

1253
    ! force unsymmetric mode when "row equilibration" is used
1254
    scaled = ListGetLogical( Solver % Values, 'Linear System Scaling', stat)
1255
    IF(.NOT.stat) scaled = .TRUE.
1256
    IF(scaled) THEN
1257
      IF(ListGetLogical( Solver % Values, 'Linear System Row Equilibration',stat)) matsym=.FALSE.
1258
    END IF
1259

1260
    IF(matsym) THEN
1261
      IF ( matspd) THEN
1262
        A % MumpsID % sym = 1
1263
      ELSE
1264
        A % MumpsID % sym = 0 ! 2=symmetric, but unsymmetric solver seems faster, at least in a few
1265
                              ! simple cases...  more testing needed...
1266
      END IF
1267
    ELSE
1268
      A % MumpsID % sym = 0
1269
    END IF
1270

1271
    CALL DMumps(A % MumpsID)
1272

1273
    IF(ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
1274

1275
    IF(ASSOCIATED(A % ParallelInfo)) THEN
1276
      n = SIZE(A % ParallelInfo % GlobalDOFs)
1277

1278
      ALLOCATE( A % Gorder(n), Owner(n) )
1279
      CALL ContinuousNumbering( A % ParallelInfo, A % Perm, A % Gorder, Owner )
1280

1281
      CALL MPI_ALLREDUCE( SUM(Owner), A % MumpsID % n, &
1282
         1, MPI_INTEGER, MPI_SUM, A % MumpsID % Comm, ierr )
1283
      DEALLOCATE(Owner)
1284
    ELSE
1285
      CALL MPI_ALLREDUCE( A % NumberOfRows, A % MumpsId % n, &
1286
            1, MPI_INTEGER, MPI_MAX, A % Comm, ierr )
1287

1288
      ALLOCATE(A % Gorder(A % NumberOFrows))
1289
      DO i=1,A % NumberOFRows
1290
        A % Gorder(i) = i
1291
      END DO
1292
    END IF
1293

1294
   ! Set matrix for Mumps (unsymmetric case)
1295
    IF (A % mumpsID % sym == 0) THEN
1296
      A % MumpsID % nz_loc = A % Rows(A % NumberOfRows+1)-1
1297

1298
      ALLOCATE( A % MumpsID % irn_loc(A % MumpsID % nz_loc) )
1299
      DO i=1,A % NumberOfRows
1300
        ip = A % Gorder(i)
1301
        DO j=A % Rows(i),A % Rows(i+1)-1
1302
          A % MumpsID % irn_loc(j) = ip
1303
        END DO
1304
      END DO
1305

1306
      ALLOCATE( A % MumpsID % jcn_loc(A % MumpsId % nz_loc) )
1307
      DO i=1,A % MumpsID % nz_loc
1308
        A % MumpsID % jcn_loc(i) = A % Gorder(A % Cols(i))
1309
      END DO
1310
      A % MumpsID % a_loc   => A % values
1311
    ELSE
1312
      ! Set matrix for Mumps (symmetric case)
1313
      nzloc = 0
1314
      DO i=1,A % NumberOfRows
1315
        ! Only output lower triangular part to Mumps
1316
        DO j=A % Rows(i),A % Diag(i)
1317
          nzloc = nzloc + 1
1318
        END DO
1319
      END DO
1320

1321
      A % MumpsID % nz_loc = nzloc
1322

1323
      ALLOCATE( A % MumpsID % irn_loc(A % MumpsID % nz_loc) )
1324
      ALLOCATE( A % MumpsID % jcn_loc(A % MumpsId % nz_loc) )
1325
      ALLOCATE( A % MumpsID % A_loc(A % MumpsId % nz_loc) )
1326

1327
      nzloc = 0
1328
      DO i=1,A % NumberOfRows
1329
        ! Only output lower triangular part to Mumps
1330
        DO j=A % Rows(i),A % Diag(i)
1331
          nzloc = nzloc + 1
1332
          A % mumpsID % IRN_loc(nzloc) = A % Gorder(i)
1333
          A % mumpsID % A_loc(nzloc) = A % Values(j)
1334
          A % mumpsID % JCN_loc(nzloc) = A % Gorder(A % Cols(j))
1335
        END DO
1336
      END DO
1337
    END IF
1338

1339

1340
    ALLOCATE(A % MumpsID % rhs(A % MumpsId % n))
1341

1342
    ! Tune verbosity of MUMPS.
1343
    i = 0
1344
    IF(InfoActive(20)) i = 1
1345
    A % MumpsID % icntl(2) = i ! suppress printing of diagnostics and warnings
1346
    A % MumpsID % icntl(3) = i ! suppress statistics
1347

1348
    A % MumpsID % icntl(4) = 1 ! the same as the two above, but doesn't seem to work.
1349
    A % MumpsID % icntl(5) = 0 ! matrix format 'assembled'
1350

1351
    icntlft = ListGetInteger(Solver % Values, &
1352
          'mumps percentage increase working space', stat)
1353
    IF (stat) THEN
1354
       A % MumpsID % icntl(14) = icntlft
1355
    END IF
1356
    A % MumpsID % icntl(18) = 3 ! 'distributed' matrix 
1357
    A % MumpsID % icntl(21) = 1 ! 'distributed' solution phase
1358

1359
    A % MumpsID % job = 4
1360
    CALL DMumps(A % MumpsID)
1361
    CALL Flush(6)
1362

1363
    A % MumpsID % lsol_loc = A % mumpsid % info(23)
1364
    ALLOCATE(A % MumpsID % sol_loc(A % MumpsId % lsol_loc))
1365
    ALLOCATE(A % MumpsID % isol_loc(A % MumpsId % lsol_loc))
1366
  END IF
1367

1368
 ! sum the rhs from all procs. Could be done for neighbours only (i guess):
1369
 ! ------------------------------------------------------------------------
1370
  A % MumpsID % RHS = 0
1371
  DO i=1,A % NumberOfRows
1372
    ip = A % Gorder(i)
1373
    A % MumpsId % RHS(ip) = b(i)
1374
  END DO
1375

1376
  ALLOCATE( dbuf(A % MumpsID % n) )
1377
  dbuf = A % MumpsId % RHS
1378
  CALL MPI_ALLREDUCE( dbuf, A % MumpsID % RHS, &
1379
    A % MumpsID % n, MPI_DOUBLE_PRECISION, MPI_SUM, A % MumpsID % Comm, ierr )
1380

1381
 ! Solution:
1382
 ! ---------
1383
  A % MumpsID % job = 3
1384
  CALL DMumps(A % MumpsID)
1385

1386
 ! Distribute the solution to all:
1387
 ! -------------------------------
1388
  A % MumpsId % Rhs = 0
1389
  DO i=1,A % MumpsID % lsol_loc
1390
    A % MumpsID % RHS(A % MumpsID % isol_loc(i)) = A % MumpsID % sol_loc(i)
1391
  END DO
1392
  dbuf = A % MumpsId % RHS
1393
  CALL MPI_ALLREDUCE( dbuf, A % MumpsID % RHS, &
1394
    A % MumpsID % N, MPI_DOUBLE_PRECISION, MPI_SUM,A %  MumpsID % Comm, ierr )
1395

1396
  DEALLOCATE(dbuf)
1397

1398
 ! Select the values which belong to us:
1399
 ! -------------------------------------
1400
  DO i=1,A % NumberOfRows
1401
    ip = A % Gorder(i)
1402
    x(i) = A % MumpsId % RHS(ip)
1403
  END DO
1404

1405
  FreeFactorize = ListGetLogical( Solver % Values, 'Linear System Free Factorization', stat )
1406
  IF ( .NOT. stat ) FreeFactorize = .TRUE.
1407
  IF ( Factorize .AND. FreeFactorize ) CALL FreeMumpsFactorizations(A)
1408

1409
#else
1410
   CALL Fatal( 'Mumps_SolveSystem', 'MUMPS Solver has not been installed.' )
1411
#endif
1412
!------------------------------------------------------------------------------
1413
  END SUBROUTINE Mumps_SolveSystem
1414
!------------------------------------------------------------------------------
1415

1416
!------------------------------------------------------------------------------
1417
!> Solves a complex linear system using MUMPS direct solver. This is a legacy solver
1418
!> with complicated dependencies. This is only available in parallel. 
1419
!------------------------------------------------------------------------------
1420
  SUBROUTINE ZMumps_SolveSystem( Solver,A,x,b )
1421
!------------------------------------------------------------------------------
1422
#ifdef HAVE_MUMPS
1423
#  if defined(ELMER_HAVE_MPI_MODULE)
1424
  USE mpi
1425
#  endif
1426
#endif
1427

1428
  TYPE(Matrix_t) :: A
1429
  TYPE(Solver_t) :: Solver
1430
  REAL(KIND=dp), TARGET :: x(*), b(*)
1431

1432
#ifdef HAVE_MUMPS
1433
#  if defined(ELMER_HAVE_MPIF_HEADER)
1434
  INCLUDE 'mpif.h'
1435
#  endif
1436

1437
  INTEGER, ALLOCATABLE :: Owner(:)
1438
  INTEGER :: i,j,k,l,n,ip,ierr,icntlft,nzloc
1439
  LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled
1440

1441
  INTEGER, ALLOCATABLE :: memb(:)
1442
  INTEGER :: Comm_active, Group_active, Group_world
1443

1444
  COMPLEX(KIND=dp), ALLOCATABLE :: dbuf(:)
1445

1446
  Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat )
1447
  IF ( .NOT. stat ) Factorize = .TRUE.
1448

1449
  IF ( Factorize .OR. .NOT.ASSOCIATED(A % ZMumpsID) ) THEN
1450
    CALL FreeMumpsFactorizations(A)
1451
    ALLOCATE(A % ZMumpsID)
1452

1453
    A % ZMumpsID % Comm = A % Comm
1454
    A % ZMumpsID % par  =  1
1455
    A % ZMumpsID % job  = -1
1456
    A % ZMumpsID % Keep =  0
1457

1458
!   matsym = ListGetLogical( Solver % Values, 'Linear System Symmetric', stat)
1459
!   matspd = ListGetLogical( Solver % Values, 'Linear System Positive Definite', stat)
1460
    matsym = .FALSE.
1461
    matspd = .FALSE.
1462

1463
    ! force unsymmetric mode when "row equilibration" is used
1464
    scaled = ListGetLogical( Solver % Values, 'Linear System Scaling', stat)
1465
    IF(.NOT.stat) scaled = .TRUE.
1466
    IF(scaled) THEN
1467
      IF(ListGetLogical( Solver % Values, 'Linear System Row Equilibration',stat)) matsym=.FALSE.
1468
    END IF
1469

1470
    A % ZMumpsID % sym = 0
1471

1472
!   IF(matsym) THEN
1473
!     IF ( matspd) THEN
1474
!       A % ZMumpsID % sym = 1
1475
!     ELSE
1476
!       A % ZMumpsID % sym = 2 ! 2=symmetric, but unsymmetric solver seems faster, at least in a few
1477
!                             ! simple cases...  more testing needed...
1478
!     END IF
1479
!   ELSE
1480
!     A % ZMumpsID % sym = 0
1481
!   END IF
1482

1483
    CALL ZMumps(A % ZMumpsID)
1484

1485
    IF(ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder)
1486

1487
    IF(ASSOCIATED(A % ParallelInfo)) THEN
1488
      n = SIZE(A % ParallelInfo % GlobalDOFs)
1489

1490
      ALLOCATE( A % Gorder(n), Owner(n) )
1491
      CALL ContinuousNumbering( A % ParallelInfo, A % Perm, A % Gorder, Owner )
1492
      CALL MPI_ALLREDUCE( SUM(Owner)/2,A % ZMumpsID % n,1,MPI_INTEGER, MPI_SUM,A % ZMumpsID % Comm,ierr )
1493
      DEALLOCATE(Owner)
1494
    ELSE
1495
      CALL MPI_ALLREDUCE( A % NumberOfRows/2, A % ZMumpsId % n, 1, MPI_INTEGER, MPI_MAX, A % Comm, ierr )
1496

1497
      ALLOCATE(A % Gorder(A % NumberOFrows))
1498
      DO i=1,A % NumberOFRows
1499
        A % Gorder(i) = i
1500
      END DO
1501
    END IF
1502

1503
   ! Set matrix for Mumps (unsymmetric case)
1504
   IF (A % ZmumpsID % sym == 0) THEN ! SYMMETRIC CASE NOT  DONE
1505
      A % ZMumpsID % nz_loc = (A % Rows(A % NumberOfRows+1)-1)/4
1506

1507
      ALLOCATE( A % ZMumpsID % irn_loc(A % ZMumpsID % nz_loc) )
1508
      ALLOCATE( A % ZMumpsID % jcn_loc(A % ZMumpsId % nz_loc) )
1509
      ALLOCATE( A % ZMumpsID % a_loc(A % ZMumpsId % nz_loc) )
1510

1511
      nzloc = 0
1512
      DO i=1,A % NumberOfRows,2
1513
        ip = (A % Gorder(i)-1)/2 + 1
1514
        DO j=A % Rows(i), A % Rows(i+1)-1,2
1515
          nzloc = nzloc + 1
1516
          A % ZMumpsID % irn_loc(nzloc) = ip
1517
          A % ZMumpsID % jcn_loc(nzloc) = (A % Gorder(A % Cols(j))-1)/2+1
1518
          A % ZMumpsID % a_loc(nzloc)   = CMPLX( A % Values(j), -A % Values(j+1), KIND=dp)
1519
        END DO
1520
      END DO
1521
    ELSE
1522
      ! Set matrix for Mumps (symmetric case)
1523
      nzloc = 0
1524
      DO i=1,A % NumberOfRows,2
1525
        ! Only output lower triangular part to Mumps
1526
        DO j=A % Rows(i),A % Diag(i),2
1527
          nzloc = nzloc + 1
1528
        END DO
1529
      END DO
1530

1531
      A % ZMumpsID % nz_loc = nzloc
1532

1533
      ALLOCATE( A % ZMumpsID % irn_loc(A % ZMumpsID % nz_loc) )
1534
      ALLOCATE( A % ZMumpsID % jcn_loc(A % ZMumpsId % nz_loc) )
1535
      ALLOCATE( A % ZMumpsID % A_loc(A % ZMumpsId % nz_loc) )
1536

1537
      nzloc = 0
1538
      DO i=1,A % NumberOfRows,2
1539
        ! Only output lower triangular part to Mumps
1540
        DO j=A % Rows(i),A % Diag(i),2
1541
          nzloc = nzloc + 1
1542
          A % ZmumpsID % IRN_loc(nzloc) = (A % Gorder(i)-1)/2+1
1543
          A % ZmumpsID % JCN_loc(nzloc) = (A % Gorder(A % Cols(j))-1)/2+1
1544
          A % ZmumpsID % A_loc(nzloc) = CMPLX( A % Values(j), -A % Values(j+1), KIND=dp)
1545
        END DO
1546
      END DO
1547
    END IF
1548

1549
    ALLOCATE(A % ZMumpsID % rhs(A % ZMumpsId % n))
1550

1551
    ! Tune verbosity of MUMPS.
1552
    i = 0
1553
    IF(InfoActive(20)) i = 1
1554
    A % ZMumpsID % icntl(2) = i ! suppress printing of diagnostics and warnings
1555
    A % ZMumpsID % icntl(3) = i ! suppress statistics
1556

1557
    A % ZMumpsID % icntl(4) = 1 ! the same as the two above, but doesn't seem to work.
1558
    A % ZMumpsID % icntl(5) = 0 ! matrix format 'assembled'
1559

1560
    icntlft = ListGetInteger(Solver % Values, 'mumps percentage increase working space', stat)
1561
    IF (stat) THEN
1562
       A % ZMumpsID % icntl(14) = icntlft
1563
    END IF
1564
    A % ZMumpsID % icntl(18) = 3 ! 'distributed' matrix 
1565
    A % ZMumpsID % icntl(21) = 1 ! 'distributed' solution phase
1566

1567
    A % ZMumpsID % job = 4
1568
    CALL ZMumps(A % ZMumpsID)
1569
    CALL Flush(6)
1570

1571
    A % ZMumpsID % lsol_loc = A % Zmumpsid % info(23)
1572
    ALLOCATE(A % ZMumpsID % sol_loc(A % ZMumpsId % lsol_loc))
1573
    ALLOCATE(A % ZMumpsID % isol_loc(A % ZMumpsId % lsol_loc))
1574
  END IF
1575

1576
 ! sum the rhs from all procs. Could be done
1577
 ! for neighbours only (i guess):
1578
 ! ------------------------------------------
1579
  A % ZMumpsID % RHS = 0
1580
  DO i=1,A % NumberOfRows,2
1581
    ip = (A % Gorder(i)-1)/2+1
1582
    A % ZMumpsId % RHS(ip) = CMPLX(b(i), b(i+1), KIND=dp)
1583
  END DO
1584
  ALLOCATE( dbuf(A % ZMumpsID % n) )
1585
  dbuf = A % ZMumpsId % RHS
1586
  CALL MPI_ALLREDUCE( dbuf, A % ZMumpsID % RHS, &
1587
    A % ZMumpsID % n, MPI_DOUBLE_COMPLEX, MPI_SUM, A % ZMumpsID % Comm, ierr )
1588

1589
 ! Solution:
1590
 ! ---------
1591
  A % ZMumpsID % job = 3
1592
  CALL ZMumps(A % ZMumpsID)
1593

1594
 ! Distribute the solution to all:
1595
 ! -------------------------------
1596
  A % ZMumpsId % Rhs = 0
1597
  DO i=1,A % ZMumpsID % lsol_loc
1598
    A % ZMumpsID % RHS(A % ZMumpsID % isol_loc(i)) = A % ZMumpsID % sol_loc(i)
1599
  END DO
1600
  dbuf = A % ZMumpsId % RHS
1601
  CALL MPI_ALLREDUCE( dbuf, A % ZMumpsID % RHS, &
1602
    A % ZMumpsID % N, MPI_DOUBLE_COMPLEX, MPI_SUM,A %  ZMumpsID % Comm, ierr )
1603

1604
  DEALLOCATE(dbuf)
1605

1606
 ! Select the values which belong to us:
1607
 ! -------------------------------------
1608
  DO i=1,A % NumberOfRows,2
1609
    ip = (A % Gorder(i)-1)/2+1
1610
    x(i)   = REAL( A % ZMumpsId % RHS(ip) )
1611
    x(i+1) = AIMAG( A % ZMumpsId % RHS(ip) )
1612
  END DO
1613

1614
  FreeFactorize = ListGetLogical( Solver % Values, 'Linear System Free Factorization', stat )
1615
  IF ( .NOT. stat ) FreeFactorize = .TRUE.
1616
  IF ( Factorize .AND. FreeFactorize ) CALL FreeMumpsFactorizations(A)
1617
#else
1618
   CALL Fatal( 'ZMumps_SolveSystem', 'MUMPS Solver has not been installed.' )
1619
#endif
1620
!------------------------------------------------------------------------------
1621
  END SUBROUTINE ZMumps_SolveSystem
1622
!------------------------------------------------------------------------------
1623

1624

1625
!------------------------------------------------------------------------------
1626
!> Solves local linear system using MUMPS direct solver. If the solved system
1627
!> is singular, optionally one possible solution is returned.
1628
!------------------------------------------------------------------------------
1629
  SUBROUTINE MumpsLocal_SolveSystem( Solver, A, x, b, Free_Fact )
1630
!------------------------------------------------------------------------------
1631
     IMPLICIT NONE
1632

1633
     TYPE(Matrix_t) :: A
1634
     TYPE(Solver_t) :: Solver
1635
     REAL(KIND=dp), TARGET :: x(*), b(*)
1636
     LOGICAL, OPTIONAL :: Free_Fact
1637

1638
     INTEGER :: i
1639
     LOGICAL :: Factorize, FreeFactorize, stat
1640

1641
#ifdef HAVE_MUMPS
1642
     ! Free local Mumps instance if requested
1643
     IF (PRESENT(Free_Fact)) THEN
1644
       IF (Free_Fact) THEN
1645
         CALL MumpsLocal_Free(A)
1646
         RETURN
1647
       END IF
1648
     END IF
1649

1650
     ! Refactorize local matrix if needed
1651
     Factorize = ListGetLogical( Solver % Values, &
1652
         'Linear System Refactorize', stat )
1653
     IF (.NOT. stat) Factorize = .TRUE.
1654

1655
     IF (Factorize .OR. .NOT. ASSOCIATED(A % mumpsIDL)) THEN
1656
       CALL MumpsLocal_Factorize(Solver, A)
1657
     END IF
1658
       
1659
     ! Set RHS
1660
     A % mumpsIDL % NRHS = 1
1661
     A % mumpsIDL % LRHS = A % mumpsIDL % n
1662
     DO i=1,A % NumberOfRows
1663
       A % mumpsIDL % RHS(i) = b(i)
1664
     END DO
1665
     ! We could use BLAS here..
1666
     ! CALL DCOPY(A % NumberOfRows, b, 1, A % mumpsIDL % RHS, 1)
1667

1668
     ! SOLUTION PHASE
1669
     A % mumpsIDL % job = 3
1670
     CALL DMumps(A % mumpsIDL)
1671

1672
     ! TODO: If solution is not local, redistribute the solution vector here
1673

1674
     ! Set local solution
1675
     DO i=1,A % NumberOfRows
1676
       x(i) = A % mumpsIDL % RHS(i)
1677
     END DO
1678
     ! We could use BLAS here..
1679
     ! CALL DCOPY(A % NumberOfRows, A % mumpsIDL % RHS, 1, x, 1)
1680

1681
     FreeFactorize = ListGetLogical( Solver % Values, &
1682
         'Linear System Free Factorization', stat )
1683
     IF (.NOT. stat) FreeFactorize = .TRUE.
1684

1685
     IF (Factorize .AND. FreeFactorize) CALL MumpsLocal_Free(A)
1686
#else
1687
     CALL Fatal( 'MumpsLocal_SolveSystem', 'MUMPS Solver has not been installed.' )
1688
#endif
1689
!------------------------------------------------------------------------------
1690
  END SUBROUTINE MumpsLocal_SolveSystem
1691
!------------------------------------------------------------------------------
1692

1693
  !------------------------------------------------------------------------------
1694
!> Solves local linear system using MUMPS direct solver. If the solved system
1695
!> is singular, optionally one possible solution is returned.
1696
!------------------------------------------------------------------------------
1697
  SUBROUTINE ZMumpsLocal_SolveSystem( Solver, A, x, b, Free_Fact )
1698
!------------------------------------------------------------------------------
1699
     IMPLICIT NONE
1700

1701
     TYPE(Matrix_t) :: A
1702
     TYPE(Solver_t) :: Solver
1703
     REAL(KIND=dp), TARGET :: x(*), b(*)
1704
     LOGICAL, OPTIONAL :: Free_Fact
1705

1706
     INTEGER :: i,j
1707
     LOGICAL :: Factorize, FreeFactorize, stat
1708

1709
#ifdef HAVE_MUMPS
1710
     ! Free local Mumps instance if requested
1711
     IF (PRESENT(Free_Fact)) THEN
1712
       IF (Free_Fact) THEN
1713
         CALL MumpsLocal_Free(A)
1714
         RETURN
1715
       END IF
1716
     END IF
1717

1718
     ! Refactorize local matrix if needed
1719
     Factorize = ListGetLogical( Solver % Values,'Linear System Refactorize', stat )
1720
     IF (.NOT. stat) Factorize = .TRUE.
1721

1722
     IF (Factorize .OR. .NOT. ASSOCIATED(A % ZmumpsIDL)) THEN
1723
       CALL ZMumpsLocal_Factorize(Solver, A)
1724
     END IF
1725

1726
     ! Set RHS
1727
     A % ZmumpsIDL % NRHS = 1
1728
     A % ZmumpsIDL % LRHS = A % ZmumpsIDL % n
1729
     j = 0
1730
     DO i=1,A % NumberOfRows,2
1731
       j = j + 1
1732
       A % ZmumpsIDL % RHS(j) = CMPLX(b(i), b(i+1), KIND=dp)
1733
     END DO
1734
     ! We could use BLAS here..
1735
     ! CALL DCOPY(A % NumberOfRows, b, 1, A % mumpsIDL % RHS, 1)
1736

1737
     ! SOLUTION PHASE
1738
     A % ZmumpsIDL % job = 3
1739
     CALL ZMumps(A % ZmumpsIDL)
1740

1741
     ! TODO: If solution is not local, redistribute the solution vector here
1742

1743
     ! Set local solution
1744
     j = 0
1745
     DO i=1,A % NumberOfRows,2
1746
       j = j + 1
1747
       x(i) = REAL(A % ZmumpsIDL % RHS(j) )
1748
       x(i+1) = AIMAG(A % ZmumpsIDL % RHS(j) )
1749
     END DO
1750
     ! We could use BLAS here..
1751
     ! CALL DCOPY(A % NumberOfRows, A % mumpsIDL % RHS, 1, x, 1)
1752

1753
     FreeFactorize = ListGetLogical( Solver % Values, 'Linear System Free Factorization', stat )
1754
     IF (.NOT. stat) FreeFactorize = .TRUE.
1755

1756
     IF (Factorize .AND. FreeFactorize) CALL MumpsLocal_Free(A)
1757
#else
1758
     CALL Fatal( 'ZMumpsLocal_SolveSystem', 'MUMPS Solver has not been installed.' )
1759
#endif
1760
!------------------------------------------------------------------------------
1761
   END SUBROUTINE ZMumpsLocal_SolveSystem
1762
!------------------------------------------------------------------------------
1763

1764
!------------------------------------------------------------------------------
1765
!> Factorize local matrix with Mumps
1766
!------------------------------------------------------------------------------
1767
  SUBROUTINE MumpsLocal_Factorize(Solver, A)
1768
!------------------------------------------------------------------------------
1769
#ifdef HAVE_MUMPS
1770
#  if defined(ELMER_HAVE_MPI_MODULE)
1771
    USE mpi
1772
#  endif
1773
#endif
1774
    IMPLICIT NONE
1775

1776
    TYPE(Solver_t) :: Solver
1777
    TYPE(Matrix_t) :: A
1778

1779
#ifdef HAVE_MUMPS
1780
#  if defined(ELMER_HAVE_MPIF_HEADER)
1781
    INCLUDE 'mpif.h'
1782
#  endif
1783

1784
    INTEGER :: i, j, n, nz, allocstat, icntlft, ptype, nzloc
1785
    LOGICAL :: matpd, matsym, nullpiv, stat
1786

1787
    ! INTEGER :: myrank, ierr
1788
    ! CHARACTER(len=32) :: buf
1789

1790
    IF ( ASSOCIATED(A % mumpsIDL) ) THEN
1791
         CALL MumpsLocal_Free(A)
1792
    END IF
1793

1794
    ALLOCATE(A % mumpsIDL)
1795

1796
    ! INITIALIZATION PHASE
1797

1798
    ! Initialize local instance of Mumps
1799
    ! TODO (Hybridization): change this if local system needs to be solved
1800
    ! with several cores
1801
    A % mumpsIDL % COMM = MPI_COMM_SELF
1802
    A % mumpsIDL % PAR = 1 ! Host (=self) takes part in factorization
1803

1804
    ! Check if matrix is symmetric or spd
1805
    matsym = ListGetLogical(Solver % Values, &
1806
                            'Linear System Symmetric', stat)
1807
    matpd = ListGetLogical(Solver % Values, &
1808
                           'Linear System Positive Definite', stat)
1809

1810
    A % mumpsIDL % SYM = 0
1811
    IF (matsym) THEN
1812
      IF (matpd) THEN
1813
        ! Matrix is symmetric positive definite
1814
        A % mumpsIDL % SYM = 1
1815
      ELSE
1816
        ! Matrix is symmetric
1817
        A % mumpsIDL % SYM = 2
1818
     END IF
1819
    ELSE
1820
      ! Matrix is unsymmetric
1821
      A % mumpsIDL % SYM = 0
1822
    END IF
1823
    A % mumpsIDL % JOB  = -1 ! Initialize
1824
    CALL DMumps(A % mumpsIDL)
1825

1826
    ! FACTORIZE PHASE
1827

1828
    ! Set stdio parameters
1829
    A % mumpsIDL % ICNTL(1)  = 6  ! Error messages to stdout
1830
    A % mumpsIDL % ICNTL(2)  = -1 ! No diagnostic and warning messages
1831
    A % mumpsIDL % ICNTL(3)  = -1 ! No statistic messages
1832
    A % mumpsIDL % ICNTL(4)  = 1  ! Print only error messages
1833

1834
    ! Set matrix format
1835
    A % mumpsIDL % ICNTL(5)  = 0  ! Assembled matrix format
1836
    A % mumpsIDL % ICNTL(18) = 0 ! Centralized matrix
1837
    A % mumpsIDL % ICNTL(21) = 0 ! Centralized dense solution phase
1838

1839
    ! Check if solution of singular systems is ok
1840
    A % mumpsIDL % ICNTL(24) = 0
1841
    nullpiv = ListGetLogical(Solver % Values, &
1842
                            'Mumps Solve Singular', stat)
1843
    IF (nullpiv) THEN
1844
      A % mumpsIDL % ICNTL(24) = 1
1845
      A % mumpsIDL % CNTL(1) = 1D-2     ! Pivoting threshold
1846
      A % mumpsIDL % CNTL(3) = 1D-9     ! Null pivot detection threshold
1847
      A % mumpsIDL % CNTL(5) = 1D6      ! Fixation value for null pivots
1848
      A % mumpsIDL % CNTL(13) = 1       ! Do not use ScaLAPACK on the root node
1849
      ! TODO: if needed, here set CNTL(3) and CNTL(5) as parameters for
1850
      ! more accurate null pivot detection
1851
    END IF
1852

1853
    ! Set permutation strategy for Mumps
1854
    ptype = ListGetInteger(Solver % Values, &
1855
                                'Mumps Permutation Type', stat)
1856
    IF (stat) THEN
1857
      A % mumpsIDL % ICNTL(6) = ptype
1858
    END IF
1859

1860
    ! TODO: Change this if system is larger than local.
1861
    ! For larger than local systems define global->local numbering
1862
    n = A % NumberofRows
1863
    nz = A % Rows(A % NumberOfRows+1)-1
1864
    A % mumpsIDL % N  = n
1865
    A % mumpsIDL % NZ = nz
1866
    ! A % mumpsIDL % nz_loc = nz
1867

1868
    ! Allocate rows and columns for MUMPS
1869
    ALLOCATE( A % mumpsIDL % IRN(nz), &
1870
          A % mumpsIDL % JCN(nz), &
1871
          A % mumpsIDL % A(nz), STAT=allocstat)
1872
    IF (allocstat /= 0) THEN
1873
      CALL Fatal('MumpsLocal_Factorize', &
1874
            'Memory allocation for MUMPS row and column indices failed.')
1875
    END IF
1876

1877
    ! Set matrix for Mumps (unsymmetric case)
1878
    IF (A % mumpsIDL % sym == 0) THEN
1879
      DO i=1,A % NumberOfRows
1880
         DO j=A % Rows(i),A % Rows(i+1)-1
1881
            A % mumpsIDL % IRN(j) = i
1882
         END DO
1883
       END DO
1884

1885
       ! Set columns and values
1886
       DO i=1,A % mumpsIDL % nz
1887
         A % mumpsIDL % JCN(i) = A % Cols(i)
1888
       END DO
1889
       DO i=1,A % mumpsIDL % nz
1890
         A % mumpsIDL % A(i) = A % Values(i)
1891
       END DO
1892
    ELSE
1893
      ! Set matrix for Mumps (symmetric case)
1894
      nzloc = 0
1895
      DO i=1,A % NumberOfRows
1896
        DO j=A % Rows(i),A % Rows(i+1)-1
1897
          ! Only output lower triangular part to Mumps
1898
          IF (i<=A % Cols(j)) THEN
1899
            nzloc = nzloc + 1
1900
            A % mumpsIDL % IRN(nzloc) = i
1901
            A % mumpsIDL % JCN(nzloc) = A % Cols(j)
1902
            A % mumpsIDL % A(nzloc) = A % Values(j)
1903
          END IF
1904
        END DO
1905
      END DO
1906
    END IF
1907

1908
    icntlft = ListGetInteger(Solver % Values, 'mumps percentage increase working space', stat)
1909
    IF (stat) THEN
1910
       A % mumpsIDL % ICNTL(14) = icntlft
1911
    END IF
1912

1913
    A % mumpsIDL % JOB = 1 ! Perform analysis
1914
    CALL DMumps(A % mumpsIDL)
1915
    CALL Flush(6)
1916

1917
    ! Check return status
1918
    IF (A % mumpsIDL % INFO(1)<0) THEN
1919
      CALL Fatal('MumpsLocal_Factorize','Mumps analysis phase failed')
1920
    END IF
1921

1922
    A % mumpsIDL % JOB = 2 ! Perform factorization
1923
    CALL DMumps(A % mumpsIDL)
1924
    CALL Flush(6)
1925

1926
    ! Check return status
1927
    IF (A % mumpsIDL % INFO(1)<0) THEN
1928
      CALL Fatal('MumpsLocal_Factorize','Mumps factorize phase failed')
1929
    END IF
1930

1931
    ! Allocate RHS
1932
    ALLOCATE(A % mumpsIDL % RHS(A % mumpsIDL % N), STAT=allocstat)
1933
    IF (allocstat /= 0) THEN
1934
         CALL Fatal('MumpsLocal_Factorize', &
1935
                 'Memory allocation for MUMPS solution vector and RHS failed.' )
1936
    END IF
1937
#else
1938
   CALL Fatal( 'MumpsLocal_Factorize', 'MUMPS Solver has not been installed.' )
1939
#endif
1940
!------------------------------------------------------------------------------
1941
  END SUBROUTINE MumpsLocal_Factorize
1942
!------------------------------------------------------------------------------
1943

1944
!------------------------------------------------------------------------------
1945
!> Factorize local matrix with Mumps
1946
!------------------------------------------------------------------------------
1947
  SUBROUTINE ZMumpsLocal_Factorize(Solver, A)
1948
!------------------------------------------------------------------------------
1949
#ifdef HAVE_MUMPS
1950
#  if defined(ELMER_HAVE_MPI_MODULE)
1951
    USE mpi
1952
#  endif
1953
#endif
1954
    IMPLICIT NONE
1955

1956
    TYPE(Solver_t) :: Solver
1957
    TYPE(Matrix_t) :: A
1958

1959
#ifdef HAVE_MUMPS
1960
#  if defined(ELMER_HAVE_MPIF_HEADER)
1961
    INCLUDE 'mpif.h'
1962
#  endif
1963

1964
    INTEGER :: i, j, k, n, nz, allocstat, icntlft, ptype, nzloc
1965
    LOGICAL :: matpd, matsym, nullpiv, stat
1966

1967
    ! INTEGER :: myrank, ierr
1968
    ! CHARACTER(len=32) :: buf
1969

1970
    IF ( ASSOCIATED(A % ZmumpsIDL) ) CALL MumpsLocal_Free(A)
1971

1972
    ! INITIALIZATION PHASE
1973
    ALLOCATE(A % ZmumpsIDL)
1974

1975
    ! Initialize local instance of Mumps
1976
    ! TODO (Hybridization): change this if local system needs to be solved
1977
    ! with several cores
1978
    A % ZmumpsIDL % COMM = MPI_COMM_SELF
1979
    A % ZmumpsIDL % PAR = 1 ! Host (=self) takes part in factorization
1980

1981
    ! Check if matrix is symmetric or spd
1982
    matsym = ListGetLogical(Solver % Values, &
1983
        'Linear System Symmetric', stat)
1984
 
1985
    matpd = ListGetLogical(Solver % Values, &
1986
        'Linear System Positive Definite', stat)
1987
 
1988
    A % ZmumpsIDL % SYM = 0
1989
!   IF (matsym) THEN
1990
!     IF (matpd) THEN
1991
!       ! Matrix is symmetric positive definite
1992
!       A % ZmumpsIDL % SYM = 1
1993
!     ELSE
1994
!       ! Matrix is symmetric
1995
!       A % ZmumpsIDL % SYM = 2
1996
!    END IF
1997
!   ELSE
1998
!     ! Matrix is unsymmetric
1999
!     A % ZmumpsIDL % SYM = 0
2000
!   END IF
2001

2002
    A % ZmumpsIDL % JOB  = -1 ! Initialize
2003
    CALL ZMumps(A % ZmumpsIDL)
2004

2005
    ! FACTORIZE PHASE
2006

2007
    ! Set stdio parameters
2008
    A % ZmumpsIDL % ICNTL(1)  = 6  ! Error messages to stdout
2009
    A % ZmumpsIDL % ICNTL(2)  = -1 ! No diagnostic and warning messages
2010
    A % ZmumpsIDL % ICNTL(3)  = -1 ! No statistic messages
2011
    A % ZmumpsIDL % ICNTL(4)  = 1  ! Print only error messages
2012

2013
    ! Set matrix format
2014
    A % ZmumpsIDL % ICNTL(5)  = 0  ! Assembled matrix format
2015
    A % ZmumpsIDL % ICNTL(18) = 0 ! Centralized matrix
2016
    A % ZmumpsIDL % ICNTL(21) = 0 ! Centralized dense solution phase
2017

2018
    ! Check if solution of singular systems is ok
2019
    A % ZmumpsIDL % ICNTL(24) = 0
2020
    nullpiv = ListGetLogical(Solver % Values, 'Mumps Solve Singular', stat)
2021
    IF (nullpiv) THEN
2022
      A % ZmumpsIDL % ICNTL(24) = 1
2023
      A % ZmumpsIDL % CNTL(1) = 1D-2     ! Pivoting threshold
2024
      A % ZmumpsIDL % CNTL(3) = 1D-9     ! Null pivot detection threshold
2025
      A % ZmumpsIDL % CNTL(5) = 1D6      ! Fixation value for null pivots
2026
      A % ZmumpsIDL % CNTL(13) = 1       ! Do not use ScaLAPACK on the root node
2027
      ! TODO: if needed, here set CNTL(3) and CNTL(5) as parameters for
2028
      ! more accurate null pivot detection
2029
    END IF
2030

2031
    ! Set permutation strategy for Mumps
2032
    ptype = ListGetInteger(Solver % Values, 'Mumps Permutation Type', stat)
2033
    IF (stat) A % ZmumpsIDL % ICNTL(6) = ptype
2034

2035
    ! TODO: Change this if system is larger than local.
2036
    ! For larger than local systems define global->local numbering
2037
    n = A % NumberofRows / 2
2038
    nz = (A % Rows(A % NumberOfRows+1)-1) / 4
2039
    A % ZmumpsIDL % N  = n
2040
    A % ZmumpsIDL % nz = nz
2041
    ! A % mumpsIDL % nz_loc = nz
2042

2043
    ! Allocate rows and columns for MUMPS
2044
    ALLOCATE( A % ZmumpsIDL % IRN(nz), &
2045
              A % ZmumpsIDL % JCN(nz), &
2046
              A % ZmumpsIDL % A(nz), STAT=allocstat)
2047
    IF (allocstat /= 0) THEN
2048
      CALL Fatal('MumpsLocal_Factorize', &
2049
            'Memory allocation for MUMPS row and column indices failed.')
2050
    END IF
2051

2052
    ! Set matrix for Mumps (unsymmetric case)
2053
    IF (A % ZmumpsIDL % sym == 0) THEN
2054
      nzloc = 0
2055
      DO i=1,A % NumberOfRows,2
2056
         DO j=A % Rows(i),A % Rows(i+1)-1,2
2057
            nzloc = nzloc + 1
2058
            A % ZmumpsIDL % IRN(nzloc) = i/2+1
2059
         END DO
2060
       END DO
2061

2062
       ! Set columns and values
2063
       nzloc = 0
2064
       DO i=1,A % NumberOfRows,2
2065
         DO j=A % Rows(i),A % Rows(i+1)-1,2
2066
           nzloc = nzloc + 1
2067
           A % ZmumpsIDL % JCN(nzloc) = A % Cols(j)/2+1
2068
           A % ZmumpsIDL % A(nzloc) = CMPLX(A % Values(j), -A % Values(j+1), KIND=dp )
2069
         END DO
2070
       END DO
2071
    ELSE
2072
      ! Set matrix for Mumps (symmetric case)
2073
      nzloc = 0
2074
      DO i=1,A % NumberOfRows,2
2075
        DO j=A % Rows(i),A % Rows(i+1)-1,2
2076
          ! Only output lower triangular part to Mumps
2077
          IF (i<=A % Cols(j)) THEN
2078
            nzloc = nzloc + 1
2079
            A % ZmumpsIDL % IRN(nzloc) = i/2+1
2080
            A % ZmumpsIDL % JCN(nzloc) = A % Cols(j)/2+1
2081
            A % ZmumpsIDL % A(nzloc) = CMPLX(A % Values(j), -A % Values(j+1),KIND=dp)
2082
          END IF
2083
        END DO
2084
      END DO
2085
    END IF
2086

2087
    icntlft = ListGetInteger(Solver % Values, 'mumps percentage increase working space', stat)
2088
    IF (stat) THEN
2089
       A % ZmumpsIDL % ICNTL(14) = icntlft
2090
    END IF
2091

2092
    A % ZmumpsIDL % JOB = 1 ! Perform analysis
2093
    CALL ZMumps(A % ZmumpsIDL)
2094
    CALL Flush(6)
2095

2096
    ! Check return status
2097
    IF (A % ZmumpsIDL % INFO(1)<0) THEN
2098
      CALL Fatal('MumpsLocal_Factorize','Mumps analysis phase failed')
2099
    END IF
2100

2101
    A % ZmumpsIDL % JOB = 2 ! Perform factorization
2102
    CALL ZMumps(A % ZmumpsIDL)
2103
    CALL Flush(6)
2104

2105
    ! Check return status
2106
    IF (A % ZmumpsIDL % INFO(1)<0) THEN
2107
      CALL Fatal('ZMumpsLocal_Factorize','Mumps factorize phase failed')
2108
    END IF
2109

2110
    ! Allocate RHS
2111
    ALLOCATE(A % ZmumpsIDL % RHS(A % ZmumpsIDL % N), STAT=allocstat)
2112
    IF (allocstat /= 0) THEN
2113
         CALL Fatal('ZMumpsLocal_Factorize', &
2114
                 'Memory allocation for MUMPS solution vector and RHS failed.' )
2115
    END IF
2116
#else
2117
   CALL Fatal( 'ZMumpsLocal_Factorize', 'MUMPS Solver has not been installed.' )
2118
#endif
2119
!------------------------------------------------------------------------------
2120
 END SUBROUTINE ZMumpsLocal_Factorize
2121
!------------------------------------------------------------------------------
2122

2123
!------------------------------------------------------------------------------
2124
!> Solve local nullspace using MUMPS direct solver. On exit, z will be
2125
!> allocated and will hold the jth local nullspace vectors as z(j,:).
2126
!------------------------------------------------------------------------------
2127
  SUBROUTINE MumpsLocal_SolveNullSpace(Solver, A, z, nz)
2128
!------------------------------------------------------------------------------
2129
#ifdef HAVE_MUMPS
2130
#  if defined(ELMER_HAVE_MPI_MODULE)
2131
      USE mpi
2132
#  endif
2133
#endif
2134
      IMPLICIT NONE
2135

2136
      TYPE(Solver_t) :: Solver
2137
      TYPE(Matrix_t) :: A
2138
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:,:), TARGET :: z
2139
      INTEGER :: nz, nrhs
2140

2141
#ifdef HAVE_MUMPS
2142
#  if defined(ELMER_HAVE_MPIF_HEADER)
2143
      INCLUDE 'mpif.h'
2144
#  endif
2145

2146
      INTEGER :: j,k,n, allocstat
2147
      LOGICAL :: Factorize, FreeFactorize, stat
2148

2149
      REAL(KIND=dp), DIMENSION(:), POINTER :: rhs_tmp, nspace
2150

2151
      Factorize = ListGetLogical( Solver % Values,&
2152
                                  'Linear System Refactorize', stat )
2153
      IF (.NOT. stat) Factorize = .TRUE.
2154

2155
      ! Refactorize local matrix if needed
2156
      IF ( Factorize .OR. (.NOT. ASSOCIATED(A % mumpsIDL)) ) THEN
2157
            CALL MumpsLocal_Factorize(Solver, A)
2158
      END IF
2159

2160
      ! Check symmetry condition
2161
      IF (.NOT. (A % mumpsIDL % sym == 1 .OR. A % mumpsIDL % sym == 2)) THEN
2162
         CALL Fatal( 'MumpsLocal_SolveNullSpace', &
2163
           'Mumps null space computation not supported for unsymmetric matrices.')
2164
      END IF
2165

2166
      ! Get the size of the local nullspace and allocate memory
2167
      ! TODO: this is incorrect, should be number of columns
2168
      n = A % NumberOfRows
2169
      nz = A % mumpsIDL % INFOG(28)
2170

2171
      IF (nz == 0) THEN
2172
            ALLOCATE(z(nz,n))
2173
          RETURN
2174
      END IF
2175

2176
      ALLOCATE(nspace(n*nz), STAT=allocstat)
2177
      IF (allocstat /= 0) THEN
2178
         CALL Fatal( 'MumpsLocal_SolveNullSpace', &
2179
               'Storage allocation for local nullspace failed.')
2180
      END IF
2181

2182
      rhs_tmp => A % mumpsIDL % RHS
2183
      nrhs = A % mumpsIDL % NRHS
2184
      A % mumpsIDL % RHS => nspace
2185
      A % mumpsIDL % NRHS = nz
2186

2187
      ! Solve for the complete local nullspace
2188
      A % mumpsIDL % JOB = 3
2189
      A % mumpsIDL % ICNTL(25) = -1
2190
      CALL DMumps(A % mumpsIDL)
2191
      ! Check return status
2192
      IF (A % mumpsIDL % INFO(1)<0) THEN
2193
          CALL Fatal('MumpsLocal_SolveNullSpace','Mumps nullspace solution failed')
2194
      END IF
2195

2196
      ! Restore pointer to local solution vector
2197
      A % mumpsIDL % ICNTL(25) = 0
2198
      A % mumpsIDL % RHS => rhs_tmp
2199
      A % mumpsIDL % NRHS = nrhs
2200

2201
      FreeFactorize = ListGetLogical( Solver % Values, &
2202
            'Linear System Free Factorization', stat )
2203
      IF (.NOT. stat) FreeFactorize = .TRUE.
2204

2205
      IF (Factorize .AND. FreeFactorize) THEN
2206
          CALL MumpsLocal_Free(A)
2207
       END IF
2208

2209
      ALLOCATE(z(nz,n), STAT=allocstat)
2210
      IF (allocstat /= 0) THEN
2211
          CALL Fatal( 'MumpsLocal_SolveNullSpace', &
2212
                           'Storage allocation for local nullspace failed.')
2213
      END IF
2214

2215
      ! Copy computed nullspace to z and deallocate nspace
2216
      DO j=1,nz
2217
            ! z is row major, cannot use DCOPY here
2218
        DO k=1,n
2219
              z(j,k)=nspace(n*(j-1)+k)
2220
        END DO
2221
      END DO
2222
      DEALLOCATE(nspace)
2223
#else
2224
   CALL Fatal( 'MumpsLocal_SolveNullSpace', 'MUMPS Solver has not been installed.' )
2225
#endif
2226
!------------------------------------------------------------------------------
2227
  END SUBROUTINE MumpsLocal_SolveNullSpace
2228
!------------------------------------------------------------------------------
2229

2230
!------------------------------------------------------------------------------
2231
!> Free local Mumps variables and solver internal allocations
2232
!------------------------------------------------------------------------------
2233
  SUBROUTINE MumpsLocal_Free(A)
2234
!------------------------------------------------------------------------------
2235
        IMPLICIT NONE
2236

2237
        TYPE(Matrix_t) :: A
2238

2239
#ifdef HAVE_MUMPS
2240
        IF (ASSOCIATED(A % mumpsIDL)) THEN
2241
          ! Deallocate Mumps structures
2242
          DEALLOCATE( A % mumpsIDL % irn, A % mumpsIDL % jcn, &
2243
              A % mumpsIDL % a, A % mumpsIDL % rhs)
2244

2245
          ! Free Mumps internal allocations
2246
          A % mumpsIDL % job = -2
2247
          CALL DMumps(A % mumpsIDL)
2248
          DEALLOCATE(A % mumpsIDL)
2249

2250
          A % mumpsIDL => NULL()
2251
        END IF
2252
        IF (ASSOCIATED(A % ZmumpsIDL)) THEN
2253
          ! Deallocate Mumps structures
2254
          DEALLOCATE( A % ZmumpsIDL % irn, A % ZmumpsIDL % jcn, &
2255
              A % ZmumpsIDL % a, A % ZmumpsIDL % rhs)
2256

2257
          ! Free Mumps internal allocations
2258
          A % ZmumpsIDL % job = -2
2259
          CALL ZMumps(A % ZmumpsIDL)
2260
          DEALLOCATE(A % ZmumpsIDL)
2261

2262
          A % ZmumpsIDL => NULL()
2263
        END IF
2264
#else
2265
     CALL Fatal( 'MumpsLocal_Free', 'MUMPS Solver has not been installed.' )
2266
#endif
2267
!------------------------------------------------------------------------------
2268
  END SUBROUTINE MumpsLocal_Free
2269
!------------------------------------------------------------------------------
2270

2271

2272

2273
!------------------------------------------------------------------------------
2274
!> Solves a linear system using SuperLU direct solver.
2275
!------------------------------------------------------------------------------
2276
  SUBROUTINE SuperLU_SolveSystem( Solver,A,x,b,Free_Fact )
2277
!------------------------------------------------------------------------------
2278
  LOGICAL, OPTIONAL :: Free_fact
2279
  TYPE(Matrix_t) :: A
2280
  TYPE(Solver_t) :: Solver
2281
  REAL(KIND=dp), TARGET :: x(*), b(*)
2282

2283
#ifdef HAVE_SUPERLU
2284
      LOGICAL :: stat, Factorize, FreeFactorize
2285
      integer :: n, nnz, nrhs, iinfo, iopt, nprocs
2286

2287
      interface
2288
        subroutine solve_superlu( iopt, nprocs, n, nnz, nrhs, values, cols, &
2289
                     rows, b, ldb, factors, iinfo )
2290
            use types
2291
            integer :: iopt, nprocs, n, nnz, nrhs, cols(*), rows(*), ldb, iinfo
2292
            real(kind=dp) :: values(*), b(*)
2293
            integer(kind=addrint) :: factors
2294
        end subroutine solve_superlu
2295
      end interface
2296

2297
      IF ( PRESENT(Free_Fact) ) THEN
2298
        IF ( Free_Fact ) THEN
2299
          IF ( A % SuperLU_Factors/= 0 ) THEN
2300
            iopt = 3
2301
            CALL Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, &
2302
               A % Cols, A % Rows, x, n, A % SuperLU_Factors, iinfo )
2303
            A % SuperLU_Factors = 0
2304
          END IF
2305
          RETURN
2306
        END IF
2307
      END IF
2308

2309
      n = A % NumberOfRows
2310
      nrhs = 1
2311
      x(1:n) = b(1:n)
2312
      nnz = A % Rows(n+1)-1
2313

2314
      nprocs = ListGetInteger( Solver % Values, & 
2315
              'Linear System Number of Threads', stat )
2316
      IF ( .NOT. stat ) nprocs = 1
2317
!
2318
      Factorize = ListGetLogical( Solver % Values, &
2319
         'Linear System Refactorize', stat )
2320
      IF ( .NOT. stat ) Factorize = .TRUE.
2321

2322
      IF ( Factorize .OR. A % SuperLU_Factors==0 ) THEN
2323

2324
        IF ( A % SuperLU_Factors/= 0 ) THEN
2325
          iopt = 3
2326
          call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, &
2327
                     A % Rows, x, n, A % SuperLU_Factors, iinfo )
2328
          A % SuperLU_Factors=0
2329
        END IF
2330

2331
        ! First, factorize the matrix. The factors are stored in *factors* handle.
2332
        iopt = 1
2333
        call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, &
2334
              A % Rows, x, n, A % SuperLU_Factors, iinfo )
2335
 
2336
        if (iinfo .eq. 0) then
2337
           write (*,*) 'Factorization succeeded'
2338
        else
2339
           write(*,*) 'INFO from factorization = ', iinfo
2340
        endif
2341
      END IF
2342

2343
      !
2344
      ! Second, solve the system using the existing factors.
2345
      iopt = 2
2346
      call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, &
2347
           A % Rows,  x, n, A % SuperLU_Factors, iinfo )
2348
!
2349
      if (iinfo .eq. 0) then
2350
         write (*,*) 'Solve succeeded'
2351
      else
2352
         write(*,*) 'INFO from triangular solve = ', iinfo
2353
      endif
2354

2355
! Last, free the storage allocated inside SuperLU
2356
      FreeFactorize = ListGetLogical( Solver % Values, &
2357
          'Linear System Free Factorization', stat )
2358
      IF ( .NOT. stat ) FreeFactorize = .TRUE.
2359

2360
      IF ( Factorize .AND. FreeFactorize ) THEN
2361
        iopt = 3
2362
        call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, &
2363
                    A % Rows, x, n, A % SuperLU_Factors, iinfo )
2364
        A % SuperLU_Factors = 0
2365
      END IF
2366
#endif      
2367
!------------------------------------------------------------------------------
2368
  END SUBROUTINE SuperLU_SolveSystem
2369
!------------------------------------------------------------------------------
2370

2371

2372
!------------------------------------------------------------------------------
2373
!> Permon solver
2374
!------------------------------------------------------------------------------
2375
  SUBROUTINE Permon_SolveSystem( Solver,A,x,b,Free_Fact )
2376
!------------------------------------------------------------------------------
2377
#ifdef HAVE_FETI4I
2378
  USE feti4i
2379
#  if defined(ELMER_HAVE_MPI_MODULE)
2380
  USE mpi
2381
#  endif
2382
#endif
2383

2384
  LOGICAL, OPTIONAL :: Free_Fact
2385
  TYPE(Matrix_t) :: A
2386
  TYPE(Solver_t) :: Solver
2387
  REAL(KIND=dp), TARGET :: x(*), b(*)
2388

2389
#ifdef HAVE_FETI4I
2390
#  if defined(ELMER_HAVE_MPIF_HEADER)
2391
  INCLUDE 'mpif.h'
2392
#  endif
2393

2394
  INTEGER, ALLOCATABLE :: Owner(:)
2395
  INTEGER :: i,j,n,nd,ip,ierr,icntlft,nzloc
2396
  LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled
2397

2398
  INTEGER :: n_dof_partition
2399

2400
  INTEGER, ALLOCATABLE :: memb(:), DirichletInds(:), Neighbours(:), IntOptions(:)
2401
  REAL(KIND=dp), ALLOCATABLE :: DirichletVals(:), RealOptions(:)
2402
  INTEGER :: Comm_active, Group_active, Group_world
2403

2404
  REAL(KIND=dp), ALLOCATABLE :: dbuf(:)
2405

2406

2407
  n_dof_partition = A % NumberOfRows
2408

2409
!  INTERFACE
2410
!     FUNCTION Permon_InitSolve(n, gnum, nd, dinds, dvals, n_n, n_ranks) RESULT(handle) BIND(c,name='permon_initsolve') 
2411
!        USE, INTRINSIC :: ISO_C_BINDING
2412
!        TYPE(C_PTR) :: handle
2413
!        INTEGER(C_INT), VALUE :: n, nd, n_n
2414
!        REAL(C_DOUBLE) :: dvals(*)
2415
!        INTEGER(C_INT) :: gnum(*), dinds(*), n_ranks(*)
2416
!     END FUNCTION Permon_Initsolve
2417
!
2418
!     SUBROUTINE Permon_Solve( handle, x, b ) BIND(c,name='permon_solve')
2419
!        USE, INTRINSIC :: ISO_C_BINDING
2420
!        REAL(C_DOUBLE) :: x(*), b(*)
2421
!        TYPE(C_PTR), VALUE :: handle
2422
!     END SUBROUTINE Permon_solve
2423
!  END INTERFACE
2424

2425
  IF ( PRESENT(Free_Fact) ) THEN
2426
    IF ( Free_Fact ) THEN
2427
      RETURN
2428
    END IF
2429
  END IF
2430

2431
  Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat )
2432
  IF ( .NOT. stat ) Factorize = .TRUE.
2433

2434
  IF ( Factorize .OR. .NOT.C_ASSOCIATED(A % PermonSolverInstance) ) THEN
2435
    IF ( C_ASSOCIATED(A % PermonSolverInstance) ) THEN
2436
       CALL Fatal( 'Permon', 're-entry not implemented' )
2437
    END IF
2438

2439
    nd = COUNT(A % ConstrainedDOF)
2440
    ALLOCATE(DirichletInds(nd), DirichletVals(nd))
2441
    j = 0
2442
    DO i=1,A % NumberOfRows
2443
      IF(A % ConstrainedDOF(i)) THEN
2444
        j = j + 1
2445
        DirichletInds(j) = i; DirichletVals(j) = A % Dvalues(i)
2446
      END IF
2447
    END DO
2448

2449
    !TODO sequential case not working
2450
    n = 0
2451
    ALLOCATE(neighbours(Parenv % PEs))
2452
    DO i=1,ParEnv % PEs
2453
      IF( ParEnv % IsNeighbour(i) .AND. i-1/=ParEnv % myPE) THEN
2454
        n = n + 1
2455
        neighbours(n) = i-1
2456
      END IF
2457
    END DO
2458

2459
    !A % PermonSolverInstance = Permon_InitSolve( SIZE(A % ParallelInfo % GlobalDOFs), &
2460
    !     A % ParallelInfo % GlobalDOFs, nd,  DirichletInds, DirichletVals, n, neighbours )
2461

2462
    IF( n_dof_partition /=  SIZE(A % ParallelInfo % GlobalDOFs) ) THEN
2463
      CALL Fatal( 'Permon', &
2464
        'inconsistency: A % NumberOfRows /=  SIZE(A % ParallelInfo % GlobalDOFs' )
2465
    END IF
2466

2467
    ALLOCATE(IntOptions(10))
2468
    ALLOCATE(RealOptions(1))
2469

2470
    CALL FETI4ISetDefaultIntegerOptions(IntOptions)
2471
    CALL FETI4ISetDefaultRealOptions(RealOptions)
2472

2473
    CALL FETI4ICreateInstance(A % PermonSolverInstance, A % PermonMatrix, &
2474
      A % NumberOfRows, b, A % ParallelInfo % GlobalDOFs, &
2475
      n, neighbours, &
2476
      nd, DirichletInds, DirichletVals, &
2477
      IntOptions, RealOptions)
2478
  END IF
2479

2480
  !CALL Permon_Solve( A % PermonSolverInstance, x, b )
2481
  CALL FETI4ISolve(A % PermonSolverInstance, n_dof_partition, x)
2482
#else
2483
   CALL Fatal( 'Permon_SolveSystem', 'Permon Solver has not been installed.' )
2484
#endif
2485
!------------------------------------------------------------------------------
2486
  END SUBROUTINE Permon_SolveSystem
2487
!------------------------------------------------------------------------------
2488

2489

2490

2491
!------------------------------------------------------------------------------
2492
!> Solves a linear system using Pardiso direct solver (from either MKL or
2493
!> official Pardiso distribution. If possible, MKL-version is used).
2494
!------------------------------------------------------------------------------
2495
  SUBROUTINE Pardiso_SolveSystem( Solver,A,x,b,Free_fact )
2496
!------------------------------------------------------------------------------
2497
    IMPLICIT NONE
2498

2499
    TYPE(Solver_t) :: Solver
2500
    TYPE(Matrix_t) :: A
2501
    REAL(KIND=dp), TARGET :: x(*), b(*)
2502
    LOGICAL, OPTIONAL :: Free_fact
2503

2504
! MKL version of Pardiso (interface is different)
2505
#if defined(HAVE_MKL)
2506
    INTERFACE
2507
      SUBROUTINE pardiso(pt, maxfct, mnum, mtype, phase, n, &
2508
                           values, rows, cols, perm, nrhs, iparm, msglvl, b, x, ierror)
2509
        USE Types
2510
        IMPLICIT NONE
2511
        REAL(KIND=dp) :: values(*), b(*), x(*)
2512
        INTEGER(KIND=AddrInt) :: pt(*)
2513
        INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror
2514
        INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*)
2515
      END SUBROUTINE pardiso
2516

2517
      SUBROUTINE pardisoinit(pt, mtype, iparm)
2518
        USE Types
2519
        IMPLICIT NONE
2520
        INTEGER(KIND=AddrInt) :: pt(*)
2521
        INTEGER :: mtype
2522
        INTEGER :: iparm(*)
2523
      END SUBROUTINE pardisoinit
2524
    END INTERFACE
2525

2526
    INTEGER maxfct, mnum, mtype, phase, n, nrhs, ierror, msglvl
2527
    INTEGER, POINTER :: Iparm(:)
2528
    INTEGER i, j, k, nz, idum(1), nzutd
2529
    LOGICAL :: Found, matsym, matpd
2530
    REAL*8  :: ddum(1)
2531

2532
    LOGICAL :: Factorize, FreeFactorize
2533
    INTEGER :: tlen, allocstat
2534
    CHARACTER(:), ALLOCATABLE :: mat_type
2535

2536
    REAL(KIND=dp), POINTER :: values(:)
2537
    INTEGER, POINTER  :: rows(:), cols(:)
2538

2539
    ! Check if system needs to be refactorized
2540
    Factorize = ListGetLogical( Solver % Values, &
2541
                  'Linear System Refactorize', Found )
2542
    IF ( .NOT. Found ) Factorize = .TRUE.
2543

2544
    ! Set matrix type for Pardiso
2545
    mat_type = ListGetString(Solver % Values,'Linear System Matrix Type',Found)
2546
    
2547
    IF (Found) THEN
2548
      SELECT CASE(mat_type)
2549
      CASE('positive definite')
2550
        mtype = 2
2551
      CASE('symmetric indefinite')
2552
        mtype = -2
2553
      CASE('structurally symmetric')
2554
        mtype = 1
2555
      CASE('nonsymmetric', 'general')
2556
        mtype = 11
2557
      CASE DEFAULT
2558
        mtype = 11
2559
      END SELECT
2560
    ELSE
2561
      ! Check if matrix is symmetric or spd
2562
      matsym = ListGetLogical(Solver % Values, &
2563
                            'Linear System Symmetric', Found)
2564

2565
      matpd = ListGetLogical(Solver % Values, &
2566
                    'Linear System Positive Definite', Found)
2567

2568
      IF (matsym) THEN
2569
        IF (matpd) THEN
2570
          ! Matrix is symmetric positive definite
2571
           mtype = 2
2572
         ELSE
2573
          ! Matrix is structurally symmetric (can't handle indefinite systems!!!!)
2574
          mtype = 1
2575
        END IF
2576
      ELSE
2577
        ! Matrix is unsymmetric
2578
        mtype = 11
2579
      END IF
2580
    END IF
2581

2582
    ! Free factorization if requested
2583
    IF ( PRESENT(Free_Fact) ) THEN
2584
      IF ( Free_Fact ) THEN
2585
        IF(ASSOCIATED(A % PardisoId)) THEN
2586
          phase     = -1           ! release internal memory
2587
          n = A % NumberOfRows
2588
          maxfct = 1
2589
          mnum   = 1
2590
          nrhs   = 1
2591
          msglvl = 0
2592
          CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2593
              ddum, idum, idum, idum, nrhs, A % PardisoParam, msglvl, ddum, ddum, ierror)
2594
          DEALLOCATE(A % PardisoId, A % PardisoParam)
2595
          A % PardisoId => NULL()
2596
          A % PardisoParam => NULL()
2597
        END IF
2598
        RETURN
2599
      END IF
2600
    END IF
2601

2602
    ! Get number of rows and number of nonzero elements
2603
    n = A % Numberofrows
2604
    nz = A % Rows(n+1)-1
2605

2606
    ! Copy upper triangular part of symmetric positive definite system
2607
    IF ( ABS(mtype) == 2 ) THEN
2608
      nzutd = 0
2609
      DO i=1,n
2610
        nzutd = nzutd + A % Rows(i+1)-A % Diag(i)
2611
      END DO
2612
      
2613
      ALLOCATE( values(nzutd), cols(nzutd), rows(n+1), STAT=allocstat)
2614
      IF (allocstat /= 0) THEN
2615
        CALL Fatal('Pardiso_SolveSystem', &
2616
           'Memory allocation for row and column indices failed')
2617
      END IF
2618

2619
      ! Copy upper triangular part of A to Pardiso structure
2620
      Rows(1)=1
2621
      DO i=1,n
2622
        ! Set up row pointers and copy values
2623
        nzutd = A % Rows(i+1)-A % Diag(i)
2624
        Rows(i+1) = Rows(i)+nzutd
2625
        DO j=0,nzutd-1
2626
          Cols(Rows(i)+j)=A % Cols(A%Diag(i)+j)
2627
          Values(Rows(i)+j)=A % Values(A%Diag(i)+j)
2628
        END DO
2629
      END DO
2630
      
2631
    ELSE
2632
      Cols => A % Cols
2633
      Rows => A % Rows
2634
      Values => A % Values
2635
    END IF
2636

2637
    ! Set up parameters
2638
    msglvl    = 0 ! Do not write out any info
2639
    maxfct    = 1 ! Set up space for 1 matrix at most
2640
    mnum      = 1 ! Matrix to use in the solution phase (1st and only one)
2641
    nrhs      = 1 ! Use only one RHS
2642
    iparm => A % PardisoParam
2643

2644
    ! Compute factorization if necessary
2645
    IF ( Factorize .OR. .NOT.ASSOCIATED(A % PardisoID) ) THEN
2646
      ! Free factorization
2647
      IF (ASSOCIATED(A % PardisoId)) THEN
2648
        phase = -1 ! release internal memory
2649
        CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, &
2650
                    ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror)
2651
        DEALLOCATE(A % PardisoId, A % PardisoParam)
2652
        A % PardisoId => NULL()
2653
        A % PardisoParam => NULL()
2654
      END IF
2655

2656
      ! Allocate pardiso internal structures
2657
      ALLOCATE(A % PardisoId(64), A % PardisoParam(64), STAT=allocstat)
2658
      IF (allocstat /= 0) THEN
2659
        CALL Fatal('Pardiso_SolveSystem', &
2660
                     'Memory allocation for Pardiso failed')
2661
      END IF
2662
      iparm => A % PardisoParam
2663
      iparm = 0
2664
      A % PardisoId = 0
2665
 
2666
      ! Set up scaling values for solver based on matrix type
2667
      CALL pardisoinit(A % PardisoId, mtype, iparm)
2668

2669
      ! Set up rest of parameters explicitly
2670
      iparm(1)=1      ! Do not use solver default parameters
2671
      iparm(2)=2      ! Minimize fill-in with nested dissection from Metis
2672
      iparm(3)=0      ! Reserved
2673
      iparm(4)=0      ! Always compute factorization
2674
      iparm(5)=0      ! No user input permutation
2675
      iparm(6)=0      ! Write solution vector to x
2676
      iparm(8)=5      ! Number of iterative refinement steps
2677
      iparm(18)=-1    ! Report nnz(L) and nnz(U)
2678
      iparm(19)=0     ! Do not report Mflops
2679
      iparm(27)=0     ! Do not check sparse matrix representation
2680
      iparm(35)=0     ! Use Fortran style indexing
2681
      iparm(60)=0     ! Use in-core version of Pardiso
2682

2683
      ! Perform analysis
2684
      phase = 11            ! Analysis
2685
      CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2686
            values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror)
2687

2688
      IF (ierror /= 0) THEN
2689
        WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror
2690
        CALL Fatal('Pardiso_SolveSystem','Error during analysis phase')
2691
      END IF
2692

2693
      ! Perform factorization
2694
      phase = 22            ! Factorization
2695
      CALL pardiso (A % pardisoId, maxfct, mnum, mtype, phase, n, &
2696
           values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror)
2697

2698
      IF (ierror /= 0) THEN
2699
        WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror
2700
        CALL Fatal('Pardiso_SolveSystem','Error during factorization phase')
2701
      END IF
2702
    END IF ! Compute factorization
2703

2704
    ! Perform solve
2705
    phase = 33            ! Solve, iterative refinement
2706
    CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2707
                values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror)
2708

2709
    IF (ierror /= 0) THEN
2710
      WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror
2711
      CALL Fatal('Pardiso_SolveSystem','Error during solve phase')
2712
    END IF
2713

2714
    ! Release memory if needed
2715
    FreeFactorize = ListGetLogical( Solver % Values, &
2716
                       'Linear System Free Factorization', Found )
2717
    IF ( .NOT. Found ) FreeFactorize = .TRUE.
2718

2719
    IF ( Factorize .AND. FreeFactorize ) THEN
2720
      phase     = -1           ! release internal memory
2721
      CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, &
2722
            ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror)
2723

2724
      DEALLOCATE(A % PardisoId, A % PardisoParam)
2725
      A % PardisoId => NULL()
2726
      A % PardisoParam => NULL()
2727
    END IF
2728
    IF (ABS(mtype) == 2) DEALLOCATE(Values, Rows, Cols)
2729

2730
! Distribution version of Pardiso
2731
#elif defined(HAVE_PARDISO)
2732
!..   All other variables
2733

2734
      INTERFACE
2735
        SUBROUTINE pardiso(pt, maxfct, mnum, mtype, phase, n, &
2736
          values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror, dparm)
2737
          USE Types
2738
          REAL(KIND=dp) :: values(*), b(*), x(*), dparm(*)
2739
          INTEGER(KIND=AddrInt) :: pt(*)
2740
          INTEGER :: idum(*), nrhs, iparm(*), msglvl, ierror
2741
          INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*)
2742
        END SUBROUTINE pardiso
2743

2744
        SUBROUTINE pardisoinit(pt,mtype,solver,iparm,dparm,ierror)
2745
          USE Types
2746
          INTEGER :: mtype, iparm(*),ierror,solver
2747
          REAL(KIND=dp) :: dparm(*)
2748
          INTEGER(KIND=AddrInt) :: pt(*)
2749
        END SUBROUTINE pardisoinit
2750
      END INTERFACE
2751

2752
      INTEGER maxfct, mnum, mtype, phase, n, nrhs, ierror, msglvl
2753
      INTEGER, POINTER :: Iparm(:)
2754
      INTEGER i, j, k, nz, idum(1)
2755
      LOGICAL :: Found, Symm, Posdef
2756
      REAL*8  waltime1, waltime2, ddum(1), dparm(64)
2757

2758
      LOGICAL :: Factorize, FreeFactorize
2759
      INTEGER :: tlen
2760
      CHARACTER(LEN=16) :: threads
2761

2762
      REAL(KIND=dp), POINTER :: values(:)
2763
      INTEGER, POINTER  :: rows(:), cols(:)
2764

2765
      IF ( PRESENT(Free_Fact) ) THEN
2766
        IF ( Free_Fact ) THEN
2767
          IF(ASSOCIATED(A % PardisoId)) THEN
2768
            phase     = -1           ! release internal memory
2769
            n = A % NumberOfRows
2770
            mtype  = 11
2771
            maxfct = 1
2772
            mnum   = 1
2773
            nrhs   = 1
2774
            msglvl = 0
2775
            CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2776
             ddum, idum, idum, idum, nrhs, A % PardisoParam, msglvl, ddum, ddum, ierror,dparm)
2777
            DEALLOCATE(A % PardisoId, A % PardisoParam)
2778
            A % PardisoId => NULL()
2779
            A % PardisoParam => NULL()
2780
          END IF
2781
          RETURN
2782
        END IF
2783
      END IF
2784

2785
!  .. Setup Pardiso control parameters und initialize the solvers
2786
!     internal address pointers. This is only necessary for the FIRST
2787
!     call of the PARDISO solver.
2788
!
2789
      Factorize = ListGetLogical( Solver % Values, &
2790
         'Linear System Refactorize', Found )
2791
      IF ( .NOT. Found ) Factorize = .TRUE.
2792

2793
      symm = ListGetLogical( Solver % Values, &
2794
         'Linear System Symmetric', Found )
2795

2796
      posdef = ListGetLogical( Solver % Values, &
2797
         'Linear System Positive Definite', Found )
2798

2799
      mtype = 11
2800

2801
      cols => a % cols
2802
      rows => a % rows
2803
      values => a % values
2804
      n = A % Numberofrows
2805

2806
      IF ( posdef ) THEN
2807
        nz = A % rows(n+1)-1
2808
        allocate( values(nz), cols(nz), rows(n+1) )
2809
        k = 1
2810
        do i=1,n
2811
         rows(i)=k
2812
         do j=a % rows(i), a % rows(i+1)-1
2813
           if ( a % cols(j)>=i .and. a % values(j) /= 0._dp ) then
2814
             cols(k) = a % cols(j)
2815
             values(k) = a % values(j)
2816
             k = k + 1
2817
           end if
2818
         end do
2819
        end do
2820
        rows(n+1)=k
2821
        mtype     = 2
2822
      ELSE IF ( symm ) THEN
2823
        mtype = 1
2824
      END IF
2825

2826
      msglvl    = 0       ! with statistical information
2827
      maxfct    = 1
2828
      mnum      = 1
2829
      nrhs      = 1
2830
      iparm => A % PardisoParam
2831

2832

2833
      IF ( Factorize .OR. .NOT.ASSOCIATED(A % PardisoID) ) THEN
2834
        IF(ASSOCIATED(A % PardisoId)) THEN
2835
          phase     = -1           ! release internal memory
2836
          CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, &
2837
           ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror,dparm)
2838
          DEALLOCATE(A % PardisoId, A % PardisoParam)
2839
          A % PardisoId => NULL()
2840
          A % PardisoParam => NULL()
2841
        END IF
2842

2843
        ALLOCATE(A % PardisoId(64), A % PardisoParam(64))
2844
        iparm => A % PardisoParam
2845
        CALL pardisoinit(A % PardisoId, mtype, 0, iparm, dparm, ierror )
2846

2847
!..     Reordering and Symbolic Factorization, This step also allocates
2848
!       all memory that is necessary for the factorization
2849
!
2850
        phase     = 11      ! only reordering and symbolic factorization
2851
        msglvl    = 0       ! with statistical information
2852
        maxfct    = 1
2853
        mnum      = 1
2854
        nrhs      = 1
2855

2856
!  ..   Numbers of Processors ( value of OMP_NUM_THREADS )
2857
        CALL envir( 'OMP_NUM_THREADS', threads, tlen )
2858

2859
        iparm(3) = 1
2860
        IF ( tlen>0 ) &
2861
          READ(threads(1:tlen),*) iparm(3)
2862
        IF (iparm(3)<=0) Iparm(3) = 1
2863

2864
        CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2865
          values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm)
2866

2867
        IF (ierror /= 0) THEN
2868
          WRITE(*,*) 'The following ERROR was detected: ', ierror
2869
          STOP EXIT_ERROR
2870
        END IF
2871

2872
!..     Factorization.
2873
        phase     = 22  ! only factorization
2874
        CALL pardiso (A % pardisoId, maxfct, mnum, mtype, phase, n, &
2875
         values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm)
2876

2877
        IF (ierror /= 0) THEN
2878
           WRITE(*,*) 'The following ERROR was detected: ', ierror
2879
          STOP EXIT_ERROR
2880
        ENDIF
2881
      END IF
2882

2883
!..   Back substitution and iterative refinement
2884
      phase     = 33  ! only factorization
2885
      iparm(8)  = 0   ! max number of iterative refinement steps
2886
                      ! (if perturbations occur in the factorization
2887
                      ! phase, two refinement steps are taken)
2888

2889
      CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, &
2890
       values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror, dparm)
2891

2892
!..   Termination and release of memory
2893
      FreeFactorize = ListGetLogical( Solver % Values, &
2894
          'Linear System Free Factorization', Found )
2895
      IF ( .NOT. Found ) FreeFactorize = .TRUE.
2896

2897
      IF ( Factorize .AND. FreeFactorize ) THEN
2898
        phase     = -1           ! release internal memory
2899
        CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, &
2900
          ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm)
2901

2902
        DEALLOCATE(A % PardisoId, A % PardisoParam)
2903
        A % PardisoId => NULL()
2904
        A % PardisoParam => NULL()
2905
      END IF
2906

2907
      IF(posdef) DEALLOCATE( values, rows, cols )
2908
#else
2909
      CALL Fatal( 'Parsido_SolveSystem', 'Pardiso solver has not been installed.' )
2910
#endif
2911

2912
!------------------------------------------------------------------------------
2913
  END SUBROUTINE Pardiso_SolveSystem
2914
!------------------------------------------------------------------------------
2915

2916
!------------------------------------------------------------------------------
2917
!> Solves a linear system using Cluster Pardiso direct solver from MKL
2918
!------------------------------------------------------------------------------
2919
  SUBROUTINE CPardiso_SolveSystem( Solver,A,x,b,Free_fact )
2920
!------------------------------------------------------------------------------
2921
    IMPLICIT NONE
2922

2923
    TYPE(Solver_t) :: Solver
2924
    TYPE(Matrix_t) :: A
2925
    REAL(KIND=dp), TARGET :: x(*), b(*)
2926
    LOGICAL, OPTIONAL :: Free_fact
2927

2928
! Cluster Pardiso
2929
#if defined(HAVE_MKL) && defined(HAVE_CPARDISO)
2930
    INTERFACE
2931
        SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, &
2932
              values, rows, cols, perm, nrhs, iparm, msglvl, b, x, comm, ierror)
2933
            USE Types
2934
            REAL(KIND=dp) :: values(*), b(*), x(*)
2935
            INTEGER(KIND=AddrInt) :: pt(*)
2936
            INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror
2937
            INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm
2938
        END SUBROUTINE cluster_sparse_solver
2939
    END INTERFACE
2940

2941
    INTEGER :: phase, n, ierror
2942
    INTEGER, POINTER :: Iparm(:)
2943
    INTEGER i, j, k, nz, nzutd, idum(1), nl, nt
2944
    LOGICAL :: Found, matsym, matpd
2945
    REAL(kind=dp) :: ddum(1)
2946
    REAL(kind=dp), POINTER, DIMENSION(:) :: dbuf
2947

2948
    LOGICAL :: Factorize, FreeFactorize
2949
    INTEGER :: tlen, allocstat
2950

2951
    REAL(KIND=dp), POINTER CONTIG :: values(:)
2952
    INTEGER, POINTER CONTIG :: rows(:), cols(:)
2953

2954
    ! Free factorization if requested
2955
    IF ( PRESENT(Free_Fact) ) THEN
2956
        IF ( Free_Fact ) THEN
2957
            CALL CPardiso_Free(A)
2958
        END IF
2959

2960
        RETURN
2961
    END IF
2962

2963
    ! Check if system needs to be refactorized
2964
    Factorize = ListGetLogical( Solver % Values, &
2965
                                'Linear System Refactorize', Found )
2966
    IF ( .NOT. Found ) Factorize = .TRUE.
2967

2968
    ! Compute factorization if necessary
2969
    IF ( Factorize .OR. .NOT.ASSOCIATED(A % CPardisoID) ) THEN
2970
        CALL CPardiso_Factorize(Solver, A)
2971
    END IF
2972

2973
    ! Get global start and end of domain
2974
    nl = A % CPardisoId % iparm(41)
2975
    nt = A % CPardisoId % iparm(42)
2976

2977
    ! Gather RHS
2978
    A % CPardisoId % rhs = 0D0
2979
    DO i=1,A % NumberOfRows
2980
        A % CPardisoId % rhs(A % Gorder(i)-nl+1) = b(i)
2981
    END DO
2982

2983
    ! Perform solve
2984
    phase = 33      ! Solve, iterative refinement
2985
    CALL cluster_sparse_solver(A % CPardisoId % ID, &
2986
          A % CPardisoId % maxfct, A % CPardisoId % mnum, &
2987
          A % CPardisoId % mtype,  phase, A % CPardisoId % n, &
2988
          A % CPardisoID % aa, A % CPardisoID % ia, A % CPardisoID % ja, idum, &
2989
          A % CPardisoId % nrhs, A % CPardisoID % iparm, &
2990
          A % CPardisoId % msglvl, A % CPardisoId % rhs, &
2991
          A % CPardisoId % x,  A % Comm, ierror)
2992

2993
    IF (ierror /= 0) THEN
2994
        WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror
2995
        CALL Fatal('CPardiso_SolveSystem','Error during solve phase')
2996
    END IF
2997

2998
    ! Distribute solution
2999
    DO i=1,A % NumberOfRows
3000
        x(i)=A % CPardisoId % x(A % Gorder(i)-nl+1)
3001
    END DO
3002

3003
    ! Release memory if needed
3004
    FreeFactorize = ListGetLogical( Solver % Values, &
3005
                                    'Linear System Free Factorization', Found )
3006
    IF ( .NOT. Found ) FreeFactorize = .TRUE.
3007

3008
    IF ( Factorize .AND. FreeFactorize ) THEN
3009
        CALL CPardiso_Free(A)
3010
    END IF
3011

3012
#else
3013
    CALL Fatal( 'CParsido_SolveSystem', 'Cluster Pardiso solver has not been installed.' )
3014
#endif
3015
!------------------------------------------------------------------------------
3016
  END SUBROUTINE CPardiso_SolveSystem
3017
!------------------------------------------------------------------------------
3018

3019
#if defined(HAVE_MKL) && defined(HAVE_CPARDISO)
3020
  SUBROUTINE CPardiso_Factorize(Solver, A)
3021
    IMPLICIT NONE
3022
    TYPE(Solver_t) :: Solver
3023
    TYPE(Matrix_t) :: A
3024

3025
    INTERFACE
3026
        SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, &
3027
              values, rows, cols, perm, nrhs, iparm, msglvl, b, x, comm, ierror)
3028
            USE Types
3029
            REAL(KIND=dp) :: values(*), b(*), x(*)
3030
            INTEGER(KIND=AddrInt) :: pt(*)
3031
            INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror
3032
            INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm
3033
        END SUBROUTINE cluster_sparse_solver
3034
    END INTERFACE
3035

3036
    LOGICAL :: matsym, matpd, Found
3037
    INTEGER :: i, j, k, rind, lrow, rptr, rsize, lind, tind
3038
    INTEGER :: allocstat, n, nz, nzutd, nl, nt, nOwned, nhalo, ierror
3039
    INTEGER :: phase, idum(1)
3040
    REAL(kind=dp) :: ddum(1)
3041
    REAL(kind=dp), DIMENSION(:), POINTER :: aa
3042
    INTEGER, DIMENSION(:), POINTER CONTIG :: iparm, ia, ja, owner, dsize, iperm, Order
3043

3044
    INTEGER :: fid
3045
    CHARACTER(:), ALLOCATABLE :: mat_type
3046

3047
    ! Free old factorization if necessary
3048
    IF (ASSOCIATED(A % CPardisoId)) THEN
3049
        CALL CPardiso_Free(A)
3050
    END IF
3051

3052
    ! Allocate Pardiso structure
3053
    ALLOCATE(A % CPardisoId, STAT=allocstat)
3054
    IF (allocstat /= 0) THEN
3055
    CALL Fatal('CPardiso_Factorize', &
3056
                   'Memory allocation for CPardiso failed')
3057
    END IF
3058
    ! Allocate control structures
3059
    ALLOCATE(A % CPardisoId% ID(64), A % CPardisoId % IParm(64), STAT=allocstat)
3060
    IF (allocstat /= 0) THEN
3061
        CALL Fatal('CPardiso_Factorize', &
3062
                   'Memory allocation for CPardiso failed')
3063
    END IF
3064

3065
    iparm => A % CPardisoId % IParm
3066
    ! Initialize control parameters and solver Id's
3067
    DO i=1,64
3068
        iparm(i)=0
3069
    END DO
3070
    DO i=1,64
3071
        A % CPardisoId % ID(i)=0
3072
    END DO
3073

3074
    ! Set matrix type for CPardiso
3075
    mat_type = ListGetString(Solver % Values,'Linear System Matrix Type',Found)
3076
    IF (Found) THEN
3077
      SELECT CASE(mat_type)
3078
      CASE('positive definite')
3079
        A % CPardisoID % mtype = 2
3080
      CASE('symmetric indefinite')
3081
        A % CPardisoID % mtype = -2
3082
      CASE('structurally symmetric')
3083
        A % CPardisoID % mtype = 1
3084
      CASE('nonsymmetric', 'general')
3085
        A % CPardisoID % mtype = 11
3086
      CASE DEFAULT
3087
        A % CPardisoID % mtype = 11
3088
      END SELECT
3089
    ELSE
3090
      ! Check if matrix is symmetric or spd
3091
      matsym = ListGetLogical(Solver % Values, &
3092
                            'Linear System Symmetric', Found)
3093

3094
      matpd = ListGetLogical(Solver % Values, &
3095
                    'Linear System Positive Definite', Found)
3096

3097
      IF (matsym) THEN
3098
        IF (matpd) THEN
3099
          ! Matrix is symmetric positive definite
3100
          A % CPardisoID % mtype = 2
3101
        ELSE
3102
          ! Matrix is structurally symmetric
3103
          A % CPardisoID % mtype = 1
3104
        END IF
3105
      ELSE
3106
        ! Matrix is nonsymmetric
3107
        A % CPardisoID % mtype = 11
3108
      END IF
3109
    END IF
3110

3111
    ! Set up continuous numbering for the whole computation domain
3112
    n = SIZE(A % ParallelInfo % GlobalDOFs)
3113
    ALLOCATE(A % Gorder(n), Owner(n), STAT=allocstat)
3114
    IF (allocstat /= 0) THEN
3115
         CALL Fatal('CPardiso_Factorize', &
3116
                    'Memory allocation for CPardiso global numbering failed')
3117
    END IF
3118
    CALL ContinuousNumbering(A % ParallelInfo, A % Perm, A % Gorder, Owner, nOwn=nOwned)
3119

3120
    ! Compute the number of global dofs
3121
    CALL MPI_ALLREDUCE(nOwned, A % CPardisoId % n, &
3122
                       1, MPI_INTEGER, MPI_SUM, A % Comm, ierror)
3123
    DEALLOCATE(Owner)
3124

3125
    ! Find bounds of domain
3126
    nl = A % Gorder(1)
3127
    nt = A % Gorder(1)
3128
    DO i=2,n
3129
        ! NOTE: Matrix is structurally symmetric
3130
        rind = A % Gorder(i)
3131
        nl = MIN(rind, nl)
3132
        nt = MAX(rind, nt)
3133
    END DO
3134

3135
    ! Allocate temp storage for global numbering
3136
    ALLOCATE(Order(n), iperm(n), STAT=allocstat)
3137
    IF (allocstat /= 0) THEN
3138
        CALL Fatal('CPardiso_Factorize', &
3139
                    'Memory allocation for CPardiso global numbering failed')
3140
    END IF
3141

3142
    ! Sort global numbering to build matrix
3143
    Order(1:n) = A % Gorder(1:n)
3144
    DO i=1,n
3145
        iperm(i)=i
3146
    END DO
3147
    CALL SortI(n, Order, iperm)
3148

3149
    ! Allocate storage for CPardiso matrix
3150
    nhalo = (nt-nl+1)-n
3151
    nz = A % Rows(A % NumberOfRows+1)-1
3152
    IF (ABS(A % CPardisoID % mtype) == 2) THEN
3153
        nzutd = ((nz-n)/2)+1 + n
3154
        ALLOCATE(A % CPardisoID % ia(nt-nl+2), &
3155
                 A % CPardisoId % ja(nzutd+nhalo), &
3156
                 A % CPardisoId % aa(nzutd+nhalo), &
3157
                 STAT=allocstat)
3158
    ELSE
3159
        ALLOCATE(A % CPardisoID % ia(nt-nl+2), &
3160
                 A % CPardisoId % ja(nz+nhalo), &
3161
                 A % CPardisoId % aa(nz+nhalo), &
3162
                STAT=allocstat)
3163
    END IF
3164
    IF (allocstat /= 0) THEN
3165
        CALL Fatal('CPardiso_Factorize', &
3166
                   'Memory allocation for CPardiso matrix failed')
3167
    END IF
3168
    ia => A % CPardisoId % ia
3169
    ja => A % CPardisoId % ja
3170
    aa => A % CPardisoId % aa
3171

3172
    ! Build distributed CRS matrix
3173
    ia(1) = 1
3174
    lrow = 1      ! Next row to add
3175
    rptr = 1      ! Pointer to next row to add, equals ia(lrow)
3176
    lind = Order(1)-1 ! Row pointer for the first round
3177
    
3178
    ! Add rows of matrix 
3179
    DO i=1,n
3180
      ! Add empty rows until the beginning of the row to add
3181
      ! (first round adds nothing due to choice of lind)
3182
      tind = Order(i)
3183
      rsize = (tind-lind)-1
3184
      
3185
      ! Put zeroes to the diagonal
3186
      DO j=1,rsize
3187
        ia(lrow+j)=rptr+j
3188
        ja(rptr+(j-1))=lind+j
3189
        aa(rptr+(j-1))=0D0
3190
      END DO
3191
      ! Set up row pointers
3192
      rptr = rptr + rsize
3193
      lrow = lrow + rsize
3194
      
3195
      ! Add next row
3196
      rind = iperm(i)
3197
      lind = A % rows(rind)
3198
      tind = A % rows(rind+1)
3199
      IF (ABS(A % CPardisoId % mtype) == 2) THEN
3200
        ! Add only upper triangular elements for symmetric matrices
3201
        rsize = 0
3202
        DO j=lind, tind-1
3203
          IF (A % Gorder(A % Cols(j)) >= Order(i)) THEN
3204
            ja(rptr+rsize)=A % Gorder(A % Cols(j))
3205
            aa(rptr+rsize)=A % values(j)
3206
            rsize = rsize + 1
3207
          END IF
3208
        END DO
3209

3210
      ELSE
3211
        rsize = tind-lind
3212
        DO j=lind, tind-1
3213
          ja(rptr+(j-lind))=A % Gorder(A % Cols(j))
3214
          aa(rptr+(j-lind))=A % values(j)
3215
        END DO
3216
      END IF
3217
        
3218
      ! Sort column indices
3219
      CALL SortF(rsize, ja(rptr:rptr+rsize), aa(rptr:rptr+rsize))
3220
        
3221
      ! Set up row pointers
3222
      rptr = rptr + rsize
3223
      lrow = lrow + 1
3224
      ia(lrow) = rptr
3225

3226
      lind = Order(i) ! Store row index for next round
3227
    END DO
3228

3229
    ! Deallocate temp storage
3230
    DEALLOCATE(Order, iperm)
3231

3232
    ! Set up parameters
3233
    A % CPardisoId % msglvl    = 0 ! Do not write out = 0 / write out = 1 info
3234
    A % CPardisoId % maxfct    = 1 ! Set up space for 1 matrix at most
3235
    A % CPardisoId % mnum      = 1 ! Matrix to use in the solution phase (1st and only one)
3236
    A % CPardisoId % nrhs      = 1 ! Use only one RHS
3237
    ALLOCATE(A % CPardisoId % rhs(nt-nl+1), &
3238
             A % CPardisoId % x(nt-nl+1), STAT=allocstat)
3239
    IF (allocstat /= 0) THEN
3240
        CALL Fatal('CPardiso_Factorize', &
3241
                   'Memory allocation for CPardiso rhs and solution vector x failed')
3242
    END IF
3243

3244
    ! Set up parameters explicitly
3245
    iparm(1)=1          ! Do not use = 1 / use = 0 solver default parameters
3246
    iparm(2)=2          ! Minimize fill-in with OpenMP nested dissection
3247
    iparm(5)=0          ! No user input permutation
3248
    iparm(6)=0          ! Write solution vector to x
3249
    iparm(8)=0          ! Number of iterative refinement steps
3250
    IF (A % CPardisoID % mtype ==11 .OR. &
3251
        A % CPardisoID % mtype == 13) THEN
3252
      iparm(10)=13      ! Perturbation value 10^-iparm(10) in case of small pivots
3253
      iparm(11)=1       ! Use scalings from symmetric weighted matching
3254
      iparm(13)=1       ! Use permutations from nonsymmetric weighted matching
3255
    ELSE
3256
      iparm(10)=8       ! Perturbation value 10^-iparm(10) in case of small pivots
3257
      iparm(11)=0       ! Do not use scalings from symmetric weighted matching
3258
      iparm(13)=0       ! Do not use permutations from symmetric weighted matching
3259
    END IF
3260
    
3261
    iparm(21)=1         ! Do not use Bunch Kaufman pivoting
3262
    iparm(27)=0         ! Do not check sparse matrix representation
3263
    iparm(28)=0         ! Use double precision
3264
    iparm(35)=0         ! Use Fortran indexing
3265

3266
    ! CPardiso matrix input format
3267
    iparm(40) = 2       ! Distributed solution phase, distributed solution vector
3268
    iparm(41) = nl      ! Beginning of solution domain
3269
    iparm(42) = nt      ! End of solution domain
3270

3271
    ! Perform analysis
3272
    phase = 11      ! Analysis
3273
    CALL cluster_sparse_solver(A % CPardisoId % ID, &
3274
          A % CPardisoId % maxfct, A % CPardisoId % mnum, &
3275
          A % CPardisoId % mtype, phase, A % CPardisoId % n, &
3276
          aa, ia, ja, idum, A % CPardisoId % nrhs, iparm, &
3277
          A % CPardisoId % msglvl, &
3278
          ddum, ddum, A % Comm, ierror)
3279
    IF (ierror /= 0) THEN
3280
        WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror
3281
        CALL Fatal('CPardiso_SolveSystem','Error during analysis phase')
3282
    END IF
3283

3284
    ! Perform factorization
3285
    phase = 22      ! Factorization
3286
    CALL cluster_sparse_solver(A % CPardisoId % ID, &
3287
          A % CPardisoId % maxfct, A % CPardisoId % mnum, &
3288
          A % CPardisoId % mtype, phase, A % CPardisoId % n, &
3289
          aa, ia, ja, idum, A % CPardisoId % nrhs, iparm, &
3290
          A % CPardisoId % msglvl, &
3291
          ddum, ddum,  A % Comm, ierror)
3292
    IF (ierror /= 0) THEN
3293
        WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror
3294
        CALL Fatal('CPardiso_SolveSystem','Error during factorization phase')
3295
    END IF
3296
  END SUBROUTINE CPardiso_Factorize
3297

3298

3299
  SUBROUTINE CPardiso_Free(A)
3300
    IMPLICIT NONE
3301

3302
    TYPE(Matrix_t) :: A
3303
    INTERFACE
3304
        SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, &
3305
                           values, rows, cols, perm, nrhs, iparm, &
3306
                           msglvl, b, x, comm, ierror)
3307
            USE Types
3308
            REAL(KIND=dp) :: values(*), b(*), x(*)
3309
            INTEGER(KIND=AddrInt) :: pt(*)
3310
            INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror
3311
            INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm
3312
        END SUBROUTINE cluster_sparse_solver
3313
    END INTERFACE
3314

3315
    INTEGER :: ierror, phase, idum(1)
3316
    REAL(kind=dp) :: ddum(1)
3317

3318
    IF(ASSOCIATED(A % CPardisoId)) THEN
3319
        phase     = -1           ! release internal memory
3320
        CALL cluster_sparse_solver(A % CPardisoId % ID, &
3321
              A % CPardisoID % maxfct, A % CPardisoID % mnum, &
3322
              A % CPardisoID % mtype, phase, A % CPardisoID % n, &
3323
              A % CPardisoId % aa, A % CPardisoId % ia, A % CPardisoId % ja, &
3324
              idum, A % CPardisoID % nrhs, A % CPardisoID % IParm, &
3325
              A % CPardisoID % msglvl, ddum, ddum, A % Comm, ierror)
3326

3327
        ! Deallocate global ordering
3328
        DEALLOCATE(A % Gorder)
3329

3330
        ! Deallocate control structures
3331
        DEALLOCATE(A % CPardisoId % ID)
3332
        DEALLOCATE(A % CPardisoID % IParm)
3333
        ! Deallocate matrix and permutation
3334
        DEALLOCATE(A % CPardisoId % ia, A % CPardisoId % ja, &
3335
                   A % CPardisoId % aa, A % CPardisoID % rhs, &
3336
                   A % CPardisoID % x)
3337
        ! Deallocate CPardiso structure
3338
        DEALLOCATE(A % CPardisoID)
3339
    END IF
3340
  END SUBROUTINE CPardiso_Free
3341
#endif
3342

3343

3344
!------------------------------------------------------------------------------
3345
  SUBROUTINE DirectSolver( A,x,b,Solver,Free_Fact )
3346
!------------------------------------------------------------------------------
3347

3348
    TYPE(Solver_t) :: Solver
3349
    REAL(KIND=dp) :: x(*),b(*)
3350
    TYPE(Matrix_t) :: A
3351
    LOGICAL, OPTIONAL :: Free_Fact
3352
!------------------------------------------------------------------------------
3353

3354
    LOGICAL :: GotIt
3355
    CHARACTER(:), ALLOCATABLE :: Method
3356
!------------------------------------------------------------------------------
3357

3358
    IF ( PRESENT(Free_Fact) ) THEN
3359
      IF ( Free_Fact ) THEN
3360
        CALL BandSolver( A, x, b, Free_Fact )
3361
        CALL ComplexBandSolver( A, x, b, Free_Fact )
3362
#ifdef HAVE_MUMPS
3363
        CALL FreeMumpsFactorizations(A)
3364
        CALL MumpsLocal_SolveSystem( Solver, A, x, b, Free_Fact )
3365
#endif
3366
#if defined(HAVE_MKL) || defined(HAVE_PARDISO)
3367
        CALL Pardiso_SolveSystem( Solver, A, x, b, Free_Fact )
3368
#endif
3369
#if defined(HAVE_MKL) && defined(HAVE_CPARDISO)
3370
        CALL CPardiso_SolveSystem( Solver, A, x, b, Free_Fact )
3371
#endif
3372
#ifdef HAVE_SUPERLU
3373
        CALL SuperLU_SolveSystem( Solver, A, x, b, Free_Fact )
3374
#endif
3375
#ifdef HAVE_UMFPACK
3376
        CALL Umfpack_SolveSystem( Solver, A, x, b, Free_Fact )
3377
#endif
3378
#ifdef HAVE_CHOLMOD
3379
        CALL SPQR_SolveSystem( Solver, A, x, b, Free_Fact )
3380
        CALL Cholmod_SolveSystem( Solver, A, x, b, Free_Fact )
3381
#endif
3382
#ifdef HAVE_FETI4I
3383
        CALL Permon_SolveSystem( Solver, A, x, b, Free_Fact )
3384
#endif
3385
        RETURN
3386
      END IF
3387
    END IF
3388

3389
    Method=ListGetString(Solver % Values,'Linear System Direct Method',GotIt)
3390
    IF ( .NOT. GotIt ) Method = 'banded'
3391
    
3392
    
3393
    CALL Info('DirectSolver','Using direct method: '//Method,Level=9)
3394

3395
#if !defined (HAVE_UMFPACK) && defined (HAVE_MUMPS)
3396
    IF ( Method == 'umfpack' .OR. Method == 'big umfpack' ) THEN
3397
      CALL Warn( 'CheckLinearSolverOptions', 'UMFPACK solver not installed, using MUMPS instead!' )
3398
      Method = 'mumps'
3399
    END IF
3400
#endif
3401

3402
    SELECT CASE(Method)
3403
      CASE( 'banded', 'symmetric banded' )
3404
        IF ( .NOT. A % Complex ) THEN
3405
           CALL BandSolver( A, x, b )
3406
        ELSE
3407
           CALL ComplexBandSolver( A, x, b )
3408
        END IF
3409

3410
      CASE( 'umfpack', 'big umfpack' )
3411
        CALL Umfpack_SolveSystem( Solver, A, x, b )
3412

3413
      CASE( 'cholmod' )
3414
        CALL Cholmod_SolveSystem( Solver, A, x, b )
3415

3416
      CASE( 'spqr' )
3417
        CALL SPQR_SolveSystem( Solver, A, x, b )
3418

3419
      CASE( 'smumps', 'cmumps' )
3420
        IF( A % Complex ) THEN
3421
          CALL CMumps_SolveSystem( Solver, A, x, b )
3422
        ELSE
3423
          CALL SMumps_SolveSystem( Solver, A, x, b )
3424
        END IF
3425

3426
      CASE( 'mumps', 'dmumps', 'zmumps' )
3427
        IF( A % Complex ) THEN
3428
          CALL ZMumps_SolveSystem( Solver, A, x, b )
3429
        ELSE
3430
          CALL Mumps_SolveSystem( Solver, A, x, b )
3431
        END IF
3432

3433
      CASE( 'mumpslocal' )
3434
        IF( A % Complex ) THEN
3435
          CALL ZMumpsLocal_SolveSystem( Solver, A, x, b )
3436
        ELSE
3437
          CALL MumpsLocal_SolveSystem( Solver, A, x, b )
3438
        END IF
3439
          
3440
      CASE( 'superlu' )
3441
        CALL SuperLU_SolveSystem( Solver, A, x, b )
3442

3443
      CASE( 'permon' )
3444
        CALL Permon_SolveSystem( Solver, A, x, b )
3445

3446
      CASE( 'pardiso' )
3447
        CALL Pardiso_SolveSystem( Solver, A, x, b )
3448

3449
      CASE( 'cpardiso' )
3450
        CALL CPardiso_SolveSystem( Solver, A, x, b )
3451

3452
      CASE DEFAULT
3453
        CALL Fatal( 'DirectSolver', 'Unknown direct solver method.' )
3454
    END SELECT
3455

3456
    ! We should be able to trust that a direct strategy will return a converged
3457
    ! linear system.
3458
    IF( ASSOCIATED( Solver % Variable ) ) THEN
3459
      Solver % Variable % LinConverged = 1
3460
    END IF
3461
    
3462
!------------------------------------------------------------------------------
3463
  END SUBROUTINE DirectSolver
3464
!------------------------------------------------------------------------------
3465

3466
END MODULE DirectSolve
3467

3468
!> \} ElmerLib
3469

3470
Product

Resources

Company
