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: Juha.Ruokolainen@csc.fi 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 47MODULE DirectSolve 48 49 USE CRSMatrix 50 USE Lists 51 USE BandMatrix 52 USE SParIterSolve 53 USE SparIterGlobals 54 55 IMPLICIT NONE 56 57CONTAINS 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(C_LONG) :: 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(C_LONG) :: 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(C_LONG) :: 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=C_LONG) :: ln, lsys 418 INTEGER(KIND=C_LONG), 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 DO i=1,n 460 LRows(i) = Rows(i)-1 461 LCols(i) = Cols(i)-1 462 END DO 463 ln = n ! TODO: Kludge: ln is AddrInt and n is regular INTEGER 464 CALL umf4_l_def( Control ) 465 CALL umf4_l_sym( ln,ln, LRows, LCols, Values, Symbolic, Control, iInfo ) 466 ELSE 467 Rows = Rows-1 468 Cols = Cols-1 469 CALL umf4def( Control ) 470 CALL umf4sym( n,n, Rows, Cols, Values, Symbolic, Control, iInfo ) 471 END IF 472 473 IF (iinfo(1)<0) THEN 474 PRINT *, 'Error occurred in umf4sym: ', iinfo(1) 475 STOP EXIT_ERROR 476 END IF 477 478 IF ( BigMode ) THEN 479 CALL umf4_l_num(LRows, LCols, Values, Symbolic, A % UMFPack_Numeric, Control, iInfo ) 480 ELSE 481 CALL umf4num( Rows, Cols, Values, Symbolic, A % UMFPack_Numeric, Control, iInfo ) 482 END IF 483 484 IF (iinfo(1)<0) THEN 485 PRINT*, 'Error occurred in umf4num: ', iinfo(1) 486 STOP EXIT_ERROR 487 ENDIF 488 489 IF ( BigMode ) THEN 490 DEALLOCATE( LRows, LCols ) 491 CALL umf4_l_fsym( Symbolic ) 492 ELSE 493 A % Rows = A % Rows+1 494 A % Cols = A % Cols+1 495 CALL umf4fsym( Symbolic ) 496 END IF 497 END IF 498 499 IF ( BigMode ) THEN 500 lsys = 2 501 CALL umf4_l_sol( lsys, x, b, A % UMFPack_Numeric, Control, iInfo ) 502 ELSE 503 sys = 2 504 CALL umf4sol( sys, x, b, A % UMFPack_Numeric, Control, iInfo ) 505 END IF 506 507 IF (iinfo(1)<0) THEN 508 PRINT*, 'Error occurred in umf4sol: ', iinfo(1) 509 STOP EXIT_ERROR 510 END IF 511 512 FreeFactorize = ListGetLogical( Solver % Values, & 513 'Linear System Free Factorization', stat ) 514 IF ( .NOT. stat ) FreeFactorize = .TRUE. 515 516 IF ( Factorize .AND. FreeFactorize ) THEN 517 IF ( BigMode ) THEN 518 CALL umf4_l_fnum(A % UMFPack_Numeric) 519 ELSE 520 CALL umf4fnum(A % UMFPack_Numeric) 521 END IF 522 A % UMFPack_Numeric = 0 523 END IF 524#else 525 CALL Fatal( 'UMFPack_SolveSystem', 'UMFPACK Solver has not been installed.' ) 526#endif 527!------------------------------------------------------------------------------ 528 END SUBROUTINE UMFPack_SolveSystem 529!------------------------------------------------------------------------------ 530 531 532!------------------------------------------------------------------------------ 533!> Solves a linear system using Cholmod multifrontal direct solver courtesy 534!> of University of Florida. 535!------------------------------------------------------------------------------ 536 SUBROUTINE Cholmod_SolveSystem( Solver,A,x,b,Free_fact) 537!------------------------------------------------------------------------------ 538 LOGICAL, OPTIONAL :: Free_Fact 539 TYPE(Matrix_t) :: A 540 TYPE(Solver_t) :: Solver 541 REAL(KIND=dp) :: x(*), b(*) 542 543 INTERFACE 544 SUBROUTINE cholmod_ffree(chol) BIND(c,NAME="cholmod_ffree") 545 USE Types 546 INTEGER(KIND=AddrInt) :: chol 547 END SUBROUTINE cholmod_ffree 548 549 FUNCTION cholmod_ffactorize(n,rows,cols,vals) RESULT(chol) BIND(c,NAME="cholmod_ffactorize") 550 USE Types 551 INTEGER :: n, Rows(*), Cols(*) 552 REAL(KIND=dp) :: Vals(*) 553 INTEGER(KIND=dp) :: chol 554 END FUNCTION cholmod_ffactorize 555 556 SUBROUTINE cholmod_fsolve(chol, n, x,b) BIND(c,NAME="cholmod_fsolve") 557 USE Types 558 REAL(KIND=dp) :: x(*), b(*) 559 INTEGER :: n 560 INTEGER(KIND=dp) :: chol 561 END SUBROUTINE cholmod_fsolve 562 END INTERFACE 563 564 LOGICAL :: Factorize, FreeFactorize, Found 565 566 REAL(KIND=dp), POINTER CONTIG :: Vals(:) 567 INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:) 568 569#ifdef HAVE_CHOLMOD 570 IF ( PRESENT(Free_Fact) ) THEN 571 IF ( Free_Fact ) THEN 572 IF ( A % Cholmod/=0 ) THEN 573 CALL cholmod_ffree(A % cholmod) 574 A % cholmod = 0 575 END IF 576 RETURN 577 END IF 578 END IF 579 580 Factorize = ListGetLogical( Solver % Values, & 581 'Linear System Refactorize', Found ) 582 IF ( .NOT. Found ) Factorize = .TRUE. 583 584 IF ( Factorize .OR. A% cholmod==0 ) THEN 585 IF ( A % cholmod/=0 ) THEN 586 CALL cholmod_ffree(A % cholmod) 587 A % cholmod = 0 588 END IF 589 590 Rows => A % Rows 591 Cols => A % Cols 592 Vals => A % Values 593 594 Rows=Rows-1; Cols=Cols-1 ! c numbering 595 A % Cholmod=cholmod_ffactorize(A % NumberOfRows, Rows, Cols, Vals) 596 Rows=Rows+1; Cols=Cols+1 ! fortran numbering 597 END IF 598 599 CALL cholmod_fsolve(A % cholmod, A % NumberOfRows, x, b); 600 601 FreeFactorize = ListGetLogical( Solver % Values, & 602 'Linear System Free Factorization', Found ) 603 IF ( .NOT. Found ) FreeFactorize = .TRUE. 604 605 IF ( Factorize .AND. FreeFactorize ) THEN 606 CALL cholmod_ffree(A % cholmod) 607 A % cholmod = 0 608 END IF 609#else 610 CALL Fatal( 'Cholmod_SolveSystem', 'Cholmod Solver has not been installed.' ) 611#endif 612!------------------------------------------------------------------------------ 613 END SUBROUTINE Cholmod_SolveSystem 614!------------------------------------------------------------------------------ 615 616 617!------------------------------------------------------------------------------ 618!> Solves a linear system using SuiteSparseQR multifrontal direct solver courtesy 619!> of University of Florida. 620!------------------------------------------------------------------------------ 621 SUBROUTINE SPQR_SolveSystem( Solver,A,x,b,Free_fact) 622!------------------------------------------------------------------------------ 623 LOGICAL, OPTIONAL :: Free_Fact 624 TYPE(Matrix_t) :: A 625 TYPE(Solver_t) :: Solver 626 REAL(KIND=dp) :: x(*), b(*) 627 628 INTEGER :: i 629 LOGICAL :: Factorize, FreeFactorize, Found 630 631 REAL(KIND=dp), POINTER CONTIG :: Vals(:) 632 INTEGER, POINTER CONTIG :: Rows(:), Cols(:), Diag(:) 633 634 INTERFACE 635 FUNCTION spqr_ffree(chol) RESULT(stat) BIND(c,NAME="spqr_ffree") 636 USE Types 637 INTEGER :: stat 638 INTEGER(KIND=AddrInt) :: chol 639 END FUNCTION spqr_ffree 640 641 FUNCTION spqr_ffactorize(n,rows,cols,vals) RESULT(chol) BIND(c,NAME="spqr_ffactorize") 642 USE Types 643 INTEGER :: n, rows(*), cols(*) 644 REAL(KIND=dp) :: vals(*) 645 INTEGER(KIND=AddrInt) :: chol 646 END FUNCTION spqr_ffactorize 647 648 SUBROUTINE spqr_fsolve(chol, n, x,b) BIND(c,NAME="spqr_fsolve") 649 USE Types 650 REAL(KIND=dp) :: x(*), b(*) 651 INTEGER :: n 652 INTEGER(KIND=AddrInt) :: chol 653 END SUBROUTINE spqr_fsolve 654 END INTERFACE 655 656#ifdef HAVE_CHOLMOD 657 IF ( PRESENT(Free_Fact) ) THEN 658 IF ( Free_Fact ) THEN 659 IF ( A % Cholmod/=0 ) THEN 660 IF(spqr_ffree(A % cholmod)==0) A % Cholmod=0 661 END IF 662 RETURN 663 END IF 664 END IF 665 666 Factorize = ListGetLogical( Solver % Values, & 667 'Linear System Refactorize', Found ) 668 IF ( .NOT. Found ) Factorize = .TRUE. 669 670 IF ( Factorize .OR. A% cholmod==0 ) THEN 671 IF ( A % cholmod/=0 ) THEN 672 i=spqr_ffree(A % cholmod) 673 A % cholmod = 0 674 END IF 675 676 Rows => A % Rows 677 Cols => A % Cols 678 Vals => A % Values 679 680 Rows=Rows-1; Cols=Cols-1 ! c numbering 681 A % Cholmod=spqr_ffactorize(A % NumberOfRows, Rows, Cols, Vals) 682 Rows=Rows+1; Cols=Cols+1 ! fortran numbering 683 END IF 684 685 CALL spqr_fsolve(A % cholmod, A % NumberOfRows, x, b); 686 687 FreeFactorize = ListGetLogical( Solver % Values, & 688 'Linear System Free Factorization', Found ) 689 IF ( .NOT. Found ) FreeFactorize = .TRUE. 690 691 IF ( Factorize .AND. FreeFactorize ) THEN 692 i=spqr_ffree(A % cholmod) 693 A % cholmod = 0 694 END IF 695#else 696 CALL Fatal( 'SPQR_SolveSystem', 'SPQR Solver has not been installed.' ) 697#endif 698!------------------------------------------------------------------------------ 699 END SUBROUTINE SPQR_SolveSystem 700!------------------------------------------------------------------------------ 701 702 703!------------------------------------------------------------------------------ 704!> Solves a linear system using MUMPS direct solver. This is a legacy solver 705!> with complicated dependencies. This is only available in parallel. 706!------------------------------------------------------------------------------ 707 SUBROUTINE Mumps_SolveSystem( Solver,A,x,b,Free_Fact ) 708!------------------------------------------------------------------------------ 709 710 LOGICAL, OPTIONAL :: Free_Fact 711 TYPE(Matrix_t) :: A 712 TYPE(Solver_t) :: Solver 713 REAL(KIND=dp), TARGET :: x(*), b(*) 714 715#ifdef HAVE_MUMPS 716 INCLUDE 'mpif.h' 717 718 INTEGER, ALLOCATABLE :: Owner(:) 719 INTEGER :: i,j,n,ip,ierr,icntlft,nzloc 720 LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled 721 722 INTEGER, ALLOCATABLE :: memb(:) 723 INTEGER :: Comm_active, Group_active, Group_world 724 725 REAL(KIND=dp), ALLOCATABLE :: dbuf(:) 726 727 728 IF ( PRESENT(Free_Fact) ) THEN 729 IF ( Free_Fact ) THEN 730 IF ( ASSOCIATED(A % MumpsID) ) THEN 731 DEALLOCATE( A % MumpsID % irn_loc, & 732 A % MumpsID % jcn_loc, A % MumpsID % rhs, & 733 A % MumpsID % isol_loc, A % MumpsID % sol_loc, A % Gorder) 734 A % Gorder=>Null() 735 736 A % MumpsID % job = -2 737 CALL DMumps(A % MumpsID) 738 DEALLOCATE(A % MumpsID) 739 END IF 740 RETURN 741 END IF 742 END IF 743 744 Factorize = ListGetLogical( Solver % Values, & 745 'Linear System Refactorize', stat ) 746 IF ( .NOT. stat ) Factorize = .TRUE. 747 748 IF ( Factorize .OR. .NOT.ASSOCIATED(A % MumpsID) ) THEN 749 IF ( ASSOCIATED(A % MumpsID) ) THEN 750 DEALLOCATE( A % MumpsID % irn_loc, & 751 A % MumpsID % jcn_loc, A % MumpsID % Rhs, & 752 A % MumpsID % isol_loc, A % MumpsID % sol_loc, A % Gorder) 753 A % MumpsID % job = -2 754 CALL DMumps(A % MumpsID) 755 DEALLOCATE(A % MumpsID) 756 END IF 757 758 ALLOCATE(A % MumpsID) 759 760 A % MumpsID % Comm = A % Comm 761 A % MumpsID % par = 1 762 A % MumpsID % job = -1 763 A % MumpsID % Keep = 0 764 765 matsym = ListGetLogical( Solver % Values, 'Linear System Symmetric', stat) 766 matspd = ListGetLogical( Solver % Values, 'Linear System Positive Definite', stat) 767 768 ! force unsymmetric mode when "row equilibration" is used 769 scaled = ListGetLogical( Solver % Values, 'Linear System Scaling', stat) 770 IF(.NOT.stat) scaled = .TRUE. 771 IF(scaled) THEN 772 IF(ListGetLogical( Solver % Values, 'Linear System Row Equilibration',stat)) matsym=.FALSE. 773 END IF 774 775 IF(matsym) THEN 776 IF ( matspd) THEN 777 A % MumpsID % sym = 1 778 ELSE 779 A % MumpsID % sym = 0 ! 2=symmetric, but unsymmetric solver seems faster, at least in a few 780 ! simple cases... more testing needed... 781 END IF 782 ELSE 783 A % MumpsID % sym = 0 784 END IF 785 786 CALL DMumps(A % MumpsID) 787 788 IF(ASSOCIATED(A % Gorder)) DEALLOCATE(A % Gorder) 789 790 IF(ASSOCIATED(A % ParallelInfo)) THEN 791 n = SIZE(A % ParallelInfo % GlobalDOFs) 792 793 ALLOCATE( A % Gorder(n), Owner(n) ) 794 CALL ContinuousNumbering( A % ParallelInfo, & 795 A % Perm, A % Gorder, Owner ) 796 797 CALL MPI_ALLREDUCE( SUM(Owner), A % MumpsID % n, & 798 1, MPI_INTEGER, MPI_SUM, A % MumpsID % Comm, ierr ) 799 DEALLOCATE(Owner) 800 ELSE 801 CALL MPI_ALLREDUCE( A % NumberOfRows, A % MumpsId % n, & 802 1, MPI_INTEGER, MPI_MAX, A % Comm, ierr ) 803 804 ALLOCATE(A % Gorder(A % NumberOFrows)) 805 DO i=1,A % NumberOFRows 806 A % Gorder(i) = i 807 END DO 808 END IF 809 810 ! Set matrix for Mumps (unsymmetric case) 811 IF (A % mumpsID % sym == 0) THEN 812 A % MumpsID % nz_loc = A % Rows(A % NumberOfRows+1)-1 813 814 ALLOCATE( A % MumpsID % irn_loc(A % MumpsID % nz_loc) ) 815 DO i=1,A % NumberOfRows 816 ip = A % Gorder(i) 817 DO j=A % Rows(i),A % Rows(i+1)-1 818 A % MumpsID % irn_loc(j) = ip 819 END DO 820 END DO 821 822 ALLOCATE( A % MumpsID % jcn_loc(A % MumpsId % nz_loc) ) 823 DO i=1,A % MumpsID % nz_loc 824 A % MumpsID % jcn_loc(i) = A % Gorder(A % Cols(i)) 825 END DO 826 A % MumpsID % a_loc => A % values 827 ELSE 828 ! Set matrix for Mumps (symmetric case) 829 nzloc = 0 830 DO i=1,A % NumberOfRows 831 ! Only output lower triangular part to Mumps 832 DO j=A % Rows(i),A % Diag(i) 833 nzloc = nzloc + 1 834 END DO 835 END DO 836 837 A % MumpsID % nz_loc = nzloc 838 839 ALLOCATE( A % MumpsID % irn_loc(A % MumpsID % nz_loc) ) 840 ALLOCATE( A % MumpsID % jcn_loc(A % MumpsId % nz_loc) ) 841 ALLOCATE( A % MumpsID % A_loc(A % MumpsId % nz_loc) ) 842 843 nzloc = 0 844 DO i=1,A % NumberOfRows 845 ! Only output lower triangular part to Mumps 846 DO j=A % Rows(i),A % Diag(i) 847 nzloc = nzloc + 1 848 A % mumpsID % IRN_loc(nzloc) = A % Gorder(i) 849 A % mumpsID % A_loc(nzloc) = A % Values(j) 850 A % mumpsID % JCN_loc(nzloc) = A % Gorder(A % Cols(j)) 851 END DO 852 END DO 853 END IF 854 855 856 ALLOCATE(A % MumpsID % rhs(A % MumpsId % n)) 857 858 A % MumpsID % icntl(2) = 0 ! suppress printing of diagnostics and warnings 859 A % MumpsID % icntl(3) = 0 ! suppress statistics 860 A % MumpsID % icntl(4) = 1 ! the same as the two above, but doesn't seem to work. 861 A % MumpsID % icntl(5) = 0 ! matrix format 'assembled' 862 863 icntlft = ListGetInteger(Solver % Values, & 864 'mumps percentage increase working space', stat) 865 IF (stat) THEN 866 A % MumpsID % icntl(14) = icntlft 867 END IF 868 A % MumpsID % icntl(18) = 3 ! 'distributed' matrix 869 A % MumpsID % icntl(21) = 1 ! 'distributed' solution phase 870 871 A % MumpsID % job = 4 872 CALL DMumps(A % MumpsID) 873 CALL Flush(6) 874 875 A % MumpsID % lsol_loc = A % mumpsid % info(23) 876 ALLOCATE(A % MumpsID % sol_loc(A % MumpsId % lsol_loc)) 877 ALLOCATE(A % MumpsID % isol_loc(A % MumpsId % lsol_loc)) 878 END IF 879 880 ! sum the rhs from all procs. Could be done 881 ! for neighbours only (i guess): 882 ! ------------------------------------------ 883 A % MumpsID % RHS = 0 884 DO i=1,A % NumberOfRows 885 ip = A % Gorder(i) 886 A % MumpsId % RHS(ip) = b(i) 887 END DO 888 ALLOCATE( dbuf(A % MumpsID % n) ) 889 dbuf = A % MumpsId % RHS 890 CALL MPI_ALLREDUCE( dbuf, A % MumpsID % RHS, & 891 A % MumpsID % n, MPI_DOUBLE_PRECISION, MPI_SUM, A % MumpsID % Comm, ierr ) 892 893 ! Solution: 894 ! --------- 895 A % MumpsID % job = 3 896 CALL DMumps(A % MumpsID) 897 898 ! Distribute the solution to all: 899 ! ------------------------------- 900 A % MumpsId % Rhs = 0 901 DO i=1,A % MumpsID % lsol_loc 902 A % MumpsID % RHS(A % MumpsID % isol_loc(i)) = & 903 A % MumpsID % sol_loc(i) 904 END DO 905 dbuf = A % MumpsId % RHS 906 CALL MPI_ALLREDUCE( dbuf, A % MumpsID % RHS, & 907 A % MumpsID % N, MPI_DOUBLE_PRECISION, MPI_SUM,A % MumpsID % Comm, ierr ) 908 909 DEALLOCATE(dbuf) 910 911 ! Select the values which belong to us: 912 ! ------------------------------------- 913 DO i=1,A % NumberOfRows 914 ip = A % Gorder(i) 915 x(i) = A % MumpsId % RHS(ip) 916 END DO 917 918 FreeFactorize = ListGetLogical( Solver % Values, & 919 'Linear System Free Factorization', stat ) 920 IF ( .NOT. stat ) FreeFactorize = .TRUE. 921 922 IF ( Factorize .AND. FreeFactorize ) THEN 923 DEALLOCATE( A % MumpsID % irn_loc, & 924 A % MumpsID % jcn_loc, A % MumpsID % Rhs, & 925 A % MumpsID % isol_loc, A % MumpsID % sol_loc, A % Gorder) 926 927 IF ( A % MumpsID % Sym/=0 ) DEALLOCATE( A % MumpsID % A_loc ) 928 929 A % Gorder=>Null() 930 931 A % MumpsID % job = -2 932 CALL DMumps(A % MumpsID) 933 DEALLOCATE(A % MumpsID) 934 A % MumpsId => NULL() 935 END IF 936 937#else 938 CALL Fatal( 'Mumps_SolveSystem', 'MUMPS Solver has not been installed.' ) 939#endif 940!------------------------------------------------------------------------------ 941 END SUBROUTINE Mumps_SolveSystem 942!------------------------------------------------------------------------------ 943 944 945!------------------------------------------------------------------------------ 946!> Solves local linear system using MUMPS direct solver. If the solved system 947!> is singular, optionally one possible solution is returned. 948!------------------------------------------------------------------------------ 949 SUBROUTINE MumpsLocal_SolveSystem( Solver, A, x, b, Free_Fact ) 950!------------------------------------------------------------------------------ 951 IMPLICIT NONE 952 953 TYPE(Matrix_t) :: A 954 TYPE(Solver_t) :: Solver 955 REAL(KIND=dp), TARGET :: x(*), b(*) 956 LOGICAL, OPTIONAL :: Free_Fact 957 958 INTEGER :: i 959 LOGICAL :: Factorize, FreeFactorize, stat 960 961#ifdef HAVE_MUMPS 962 ! Free local Mumps instance if requested 963 IF (PRESENT(Free_Fact)) THEN 964 IF (Free_Fact) THEN 965 CALL MumpsLocal_Free(A) 966 RETURN 967 END IF 968 END IF 969 970 ! Refactorize local matrix if needed 971 Factorize = ListGetLogical( Solver % Values, & 972 'Linear System Refactorize', stat ) 973 IF (.NOT. stat) Factorize = .TRUE. 974 975 IF (Factorize .OR. .NOT. ASSOCIATED(A % mumpsIDL)) THEN 976 CALL MumpsLocal_Factorize(Solver, A) 977 END IF 978 979 ! Set RHS 980 A % mumpsIDL % NRHS = 1 981 A % mumpsIDL % LRHS = A % mumpsIDL % n 982 DO i=1,A % NumberOfRows 983 A % mumpsIDL % RHS(i) = b(i) 984 END DO 985 ! We could use BLAS here.. 986 ! CALL DCOPY(A % NumberOfRows, b, 1, A % mumpsIDL % RHS, 1) 987 988 ! SOLUTION PHASE 989 A % mumpsIDL % job = 3 990 CALL DMumps(A % mumpsIDL) 991 992 ! TODO: If solution is not local, redistribute the solution vector here 993 994 ! Set local solution 995 DO i=1,A % NumberOfRows 996 x(i)=A % mumpsIDL % RHS(i) 997 END DO 998 ! We could use BLAS here.. 999 ! CALL DCOPY(A % NumberOfRows, A % mumpsIDL % RHS, 1, x, 1) 1000 1001 FreeFactorize = ListGetLogical( Solver % Values, & 1002 'Linear System Free Factorization', stat ) 1003 IF (.NOT. stat) FreeFactorize = .TRUE. 1004 1005 IF (Factorize .AND. FreeFactorize) CALL MumpsLocal_Free(A) 1006#else 1007 CALL Fatal( 'MumpsLocal_SolveSystem', 'MUMPS Solver has not been installed.' ) 1008#endif 1009!------------------------------------------------------------------------------ 1010 END SUBROUTINE MumpsLocal_SolveSystem 1011!------------------------------------------------------------------------------ 1012 1013!------------------------------------------------------------------------------ 1014!> Factorize local matrix with Mumps 1015!------------------------------------------------------------------------------ 1016 SUBROUTINE MumpsLocal_Factorize(Solver, A) 1017!------------------------------------------------------------------------------ 1018 IMPLICIT NONE 1019 1020 TYPE(Solver_t) :: Solver 1021 TYPE(Matrix_t) :: A 1022 1023#ifdef HAVE_MUMPS 1024 INCLUDE 'mpif.h' 1025 1026 INTEGER :: i, j, n, nz, allocstat, icntlft, ptype, nzloc 1027 LOGICAL :: matpd, matsym, nullpiv, stat 1028 1029 ! INTEGER :: myrank, ierr 1030 ! CHARACTER(len=32) :: buf 1031 1032 IF ( ASSOCIATED(A % mumpsIDL) ) THEN 1033 CALL MumpsLocal_Free(A) 1034 END IF 1035 1036 ALLOCATE(A % mumpsIDL) 1037 1038 ! INITIALIZATION PHASE 1039 1040 ! Initialize local instance of Mumps 1041 ! TODO (Hybridization): change this if local system needs to be solved 1042 ! with several cores 1043 A % mumpsIDL % COMM = MPI_COMM_SELF 1044 A % mumpsIDL % PAR = 1 ! Host (=self) takes part in factorization 1045 1046 ! Check if matrix is symmetric or spd 1047 matsym = ListGetLogical(Solver % Values, & 1048 'Linear System Symmetric', stat) 1049 matpd = ListGetLogical(Solver % Values, & 1050 'Linear System Positive Definite', stat) 1051 1052 A % mumpsIDL % SYM = 0 1053 IF (matsym) THEN 1054 IF (matpd) THEN 1055 ! Matrix is symmetric positive definite 1056 A % mumpsIDL % SYM = 1 1057 ELSE 1058 ! Matrix is symmetric 1059 A % mumpsIDL % SYM = 2 1060 END IF 1061 ELSE 1062 ! Matrix is unsymmetric 1063 A % mumpsIDL % SYM = 0 1064 END IF 1065 A % mumpsIDL % JOB = -1 ! Initialize 1066 CALL DMumps(A % mumpsIDL) 1067 1068 ! FACTORIZE PHASE 1069 1070 ! Set stdio parameters 1071 A % mumpsIDL % ICNTL(1) = 6 ! Error messages to stdout 1072 A % mumpsIDL % ICNTL(2) = -1 ! No diagnostic and warning messages 1073 A % mumpsIDL % ICNTL(3) = -1 ! No statistic messages 1074 A % mumpsIDL % ICNTL(4) = 1 ! Print only error messages 1075 1076 ! Set matrix format 1077 A % mumpsIDL % ICNTL(5) = 0 ! Assembled matrix format 1078 A % mumpsIDL % ICNTL(18) = 0 ! Centralized matrix 1079 A % mumpsIDL % ICNTL(21) = 0 ! Centralized dense solution phase 1080 1081 ! Check if solution of singular systems is ok 1082 A % mumpsIDL % ICNTL(24) = 0 1083 nullpiv = ListGetLogical(Solver % Values, & 1084 'Mumps Solve Singular', stat) 1085 IF (nullpiv) THEN 1086 A % mumpsIDL % ICNTL(24) = 1 1087 A % mumpsIDL % CNTL(1) = 1D-2 ! Pivoting threshold 1088 A % mumpsIDL % CNTL(3) = 1D-9 ! Null pivot detection threshold 1089 A % mumpsIDL % CNTL(5) = 1D6 ! Fixation value for null pivots 1090 A % mumpsIDL % CNTL(13) = 1 ! Do not use ScaLAPACK on the root node 1091 ! TODO: if needed, here set CNTL(3) and CNTL(5) as parameters for 1092 ! more accurate null pivot detection 1093 END IF 1094 1095 ! Set permutation strategy for Mumps 1096 ptype = ListGetInteger(Solver % Values, & 1097 'Mumps Permutation Type', stat) 1098 IF (stat) THEN 1099 A % mumpsIDL % ICNTL(6) = ptype 1100 END IF 1101 1102 ! TODO: Change this if system is larger than local. 1103 ! For larger than local systems define global->local numbering 1104 n = A % NumberofRows 1105 nz = A % Rows(A % NumberOfRows+1)-1 1106 A % mumpsIDL % N = n 1107 A % mumpsIDL % NZ = nz 1108 ! A % mumpsIDL % nz_loc = nz 1109 1110 ! Allocate rows and columns for MUMPS 1111 ALLOCATE( A % mumpsIDL % IRN(nz), & 1112 A % mumpsIDL % JCN(nz), & 1113 A % mumpsIDL % A(nz), STAT=allocstat) 1114 IF (allocstat /= 0) THEN 1115 CALL Fatal('MumpsLocal_Factorize', & 1116 'Memory allocation for MUMPS row and column indices failed.') 1117 END IF 1118 1119 ! Set matrix for Mumps (unsymmetric case) 1120 IF (A % mumpsIDL % sym == 0) THEN 1121 DO i=1,A % NumberOfRows 1122 DO j=A % Rows(i),A % Rows(i+1)-1 1123 A % mumpsIDL % IRN(j) = i 1124 END DO 1125 END DO 1126 1127 ! Set columns and values 1128 DO i=1,A % mumpsIDL % nz 1129 A % mumpsIDL % JCN(i) = A % Cols(i) 1130 END DO 1131 DO i=1,A % mumpsIDL % nz 1132 A % mumpsIDL % A(i) = A % Values(i) 1133 END DO 1134 ELSE 1135 ! Set matrix for Mumps (symmetric case) 1136 nzloc = 0 1137 DO i=1,A % NumberOfRows 1138 DO j=A % Rows(i),A % Rows(i+1)-1 1139 ! Only output lower triangular part to Mumps 1140 IF (i<=A % Cols(j)) THEN 1141 nzloc = nzloc + 1 1142 A % mumpsIDL % IRN(nzloc) = i 1143 A % mumpsIDL % JCN(nzloc) = A % Cols(j) 1144 A % mumpsIDL % A(nzloc) = A % Values(j) 1145 END IF 1146 END DO 1147 END DO 1148 END IF 1149 1150 icntlft = ListGetInteger(Solver % Values, 'mumps percentage increase working space', stat) 1151 IF (stat) THEN 1152 A % mumpsIDL % ICNTL(14) = icntlft 1153 END IF 1154 1155 A % mumpsIDL % JOB = 1 ! Perform analysis 1156 CALL DMumps(A % mumpsIDL) 1157 CALL Flush(6) 1158 1159 ! Check return status 1160 IF (A % mumpsIDL % INFO(1)<0) THEN 1161 CALL Fatal('MumpsLocal_Factorize','Mumps analysis phase failed') 1162 END IF 1163 1164 A % mumpsIDL % JOB = 2 ! Perform factorization 1165 CALL DMumps(A % mumpsIDL) 1166 CALL Flush(6) 1167 1168 ! Check return status 1169 IF (A % mumpsIDL % INFO(1)<0) THEN 1170 CALL Fatal('MumpsLocal_Factorize','Mumps factorize phase failed') 1171 END IF 1172 1173 ! Allocate RHS 1174 ALLOCATE(A % mumpsIDL % RHS(A % mumpsIDL % N), STAT=allocstat) 1175 IF (allocstat /= 0) THEN 1176 CALL Fatal('MumpsLocal_Factorize', & 1177 'Memory allocation for MUMPS solution vector and RHS failed.' ) 1178 END IF 1179#else 1180 CALL Fatal( 'MumpsLocal_Factorize', 'MUMPS Solver has not been installed.' ) 1181#endif 1182!------------------------------------------------------------------------------ 1183 END SUBROUTINE MumpsLocal_Factorize 1184!------------------------------------------------------------------------------ 1185 1186!------------------------------------------------------------------------------ 1187!> Solve local nullspace using MUMPS direct solver. On exit, z will be 1188!> allocated and will hold the jth local nullspace vectors as z(j,:). 1189!------------------------------------------------------------------------------ 1190 SUBROUTINE MumpsLocal_SolveNullSpace(Solver, A, z, nz) 1191!------------------------------------------------------------------------------ 1192 IMPLICIT NONE 1193 1194 TYPE(Solver_t) :: Solver 1195 TYPE(Matrix_t) :: A 1196 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:,:), TARGET :: z 1197 INTEGER :: nz, nrhs 1198 1199#ifdef HAVE_MUMPS 1200 INCLUDE 'mpif.h' 1201 1202 INTEGER :: j,k,n, allocstat 1203 LOGICAL :: Factorize, FreeFactorize, stat 1204 1205 REAL(KIND=dp), DIMENSION(:), POINTER :: rhs_tmp, nspace 1206 1207 Factorize = ListGetLogical( Solver % Values,& 1208 'Linear System Refactorize', stat ) 1209 IF (.NOT. stat) Factorize = .TRUE. 1210 1211 ! Refactorize local matrix if needed 1212 IF ( Factorize .OR. (.NOT. ASSOCIATED(A % mumpsIDL)) ) THEN 1213 CALL MumpsLocal_Factorize(Solver, A) 1214 END IF 1215 1216 ! Check symmetry condition 1217 IF (.NOT. (A % mumpsIDL % sym == 1 .OR. A % mumpsIDL % sym == 2)) THEN 1218 CALL Fatal( 'MumpsLocal_SolveNullSpace', & 1219 'Mumps null space computation not supported for unsymmetric matrices.') 1220 END IF 1221 1222 ! Get the size of the local nullspace and allocate memory 1223 ! TODO: this is incorrect, should be number of columns 1224 n = A % NumberOfRows 1225 nz = A % mumpsIDL % INFOG(28) 1226 1227 IF (nz == 0) THEN 1228 ALLOCATE(z(nz,n)) 1229 RETURN 1230 END IF 1231 1232 ALLOCATE(nspace(n*nz), STAT=allocstat) 1233 IF (allocstat /= 0) THEN 1234 CALL Fatal( 'MumpsLocal_SolveNullSpace', & 1235 'Storage allocation for local nullspace failed.') 1236 END IF 1237 1238 rhs_tmp => A % mumpsIDL % RHS 1239 nrhs = A % mumpsIDL % NRHS 1240 A % mumpsIDL % RHS => nspace 1241 A % mumpsIDL % NRHS = nz 1242 1243 ! Solve for the complete local nullspace 1244 A % mumpsIDL % JOB = 3 1245 A % mumpsIDL % ICNTL(25) = -1 1246 CALL DMumps(A % mumpsIDL) 1247 ! Check return status 1248 IF (A % mumpsIDL % INFO(1)<0) THEN 1249 CALL Fatal('MumpsLocal_SolveNullSpace','Mumps nullspace solution failed') 1250 END IF 1251 1252 ! Restore pointer to local solution vector 1253 A % mumpsIDL % ICNTL(25) = 0 1254 A % mumpsIDL % RHS => rhs_tmp 1255 A % mumpsIDL % NRHS = nrhs 1256 1257 FreeFactorize = ListGetLogical( Solver % Values, & 1258 'Linear System Free Factorization', stat ) 1259 IF (.NOT. stat) FreeFactorize = .TRUE. 1260 1261 IF (Factorize .AND. FreeFactorize) THEN 1262 CALL MumpsLocal_Free(A) 1263 END IF 1264 1265 ALLOCATE(z(nz,n), STAT=allocstat) 1266 IF (allocstat /= 0) THEN 1267 CALL Fatal( 'MumpsLocal_SolveNullSpace', & 1268 'Storage allocation for local nullspace failed.') 1269 END IF 1270 1271 ! Copy computed nullspace to z and deallocate nspace 1272 DO j=1,nz 1273 ! z is row major, cannot use DCOPY here 1274 DO k=1,n 1275 z(j,k)=nspace(n*(j-1)+k) 1276 END DO 1277 END DO 1278 DEALLOCATE(nspace) 1279#else 1280 CALL Fatal( 'MumpsLocal_SolveNullSpace', 'MUMPS Solver has not been installed.' ) 1281#endif 1282!------------------------------------------------------------------------------ 1283 END SUBROUTINE MumpsLocal_SolveNullSpace 1284!------------------------------------------------------------------------------ 1285 1286!------------------------------------------------------------------------------ 1287!> Free local Mumps variables and solver internal allocations 1288!------------------------------------------------------------------------------ 1289 SUBROUTINE MumpsLocal_Free(A) 1290!------------------------------------------------------------------------------ 1291 IMPLICIT NONE 1292 1293 TYPE(Matrix_t) :: A 1294 1295#ifdef HAVE_MUMPS 1296 IF (ASSOCIATED(A % mumpsIDL)) THEN 1297 ! Deallocate Mumps structures 1298 DEALLOCATE( A % mumpsIDL % irn, A % mumpsIDL % jcn, & 1299 A % mumpsIDL % a, A % mumpsIDL % rhs) 1300 1301 ! Free Mumps internal allocations 1302 A % mumpsIDL % job = -2 1303 CALL DMumps(A % mumpsIDL) 1304 DEALLOCATE(A % mumpsIDL) 1305 1306 A % mumpsIDL => Null() 1307 END IF 1308#else 1309 CALL Fatal( 'MumpsLocal_Free', 'MUMPS Solver has not been installed.' ) 1310#endif 1311!------------------------------------------------------------------------------ 1312 END SUBROUTINE MumpsLocal_Free 1313!------------------------------------------------------------------------------ 1314 1315 1316 1317!------------------------------------------------------------------------------ 1318!> Solves a linear system using SuperLU direct solver. 1319!------------------------------------------------------------------------------ 1320 SUBROUTINE SuperLU_SolveSystem( Solver,A,x,b,Free_Fact ) 1321!------------------------------------------------------------------------------ 1322 LOGICAL, OPTIONAL :: Free_fact 1323 TYPE(Matrix_t) :: A 1324 TYPE(Solver_t) :: Solver 1325 REAL(KIND=dp), TARGET :: x(*), b(*) 1326 1327#ifdef HAVE_SUPERLU 1328 LOGICAL :: stat, Factorize, FreeFactorize 1329 integer :: n, nnz, nrhs, iinfo, iopt, nprocs 1330 1331 interface 1332 subroutine solve_superlu( iopt, nprocs, n, nnz, nrhs, values, cols, & 1333 rows, b, ldb, factors, iinfo ) 1334 use types 1335 integer :: iopt, nprocs, n, nnz, nrhs, cols(*), rows(*), ldb, iinfo 1336 real(kind=dp) :: values(*), b(*) 1337 integer(kind=addrint) :: factors 1338 end subroutine solve_superlu 1339 end interface 1340 1341 IF ( PRESENT(Free_Fact) ) THEN 1342 IF ( Free_Fact ) THEN 1343 IF ( A % SuperLU_Factors/= 0 ) THEN 1344 iopt = 3 1345 CALL Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, & 1346 A % Cols, A % Rows, x, n, A % SuperLU_Factors, iinfo ) 1347 A % SuperLU_Factors = 0 1348 END IF 1349 RETURN 1350 END IF 1351 END IF 1352 1353 n = A % NumberOfRows 1354 nrhs = 1 1355 x(1:n) = b(1:n) 1356 nnz = A % Rows(n+1)-1 1357 1358 nprocs = ListGetInteger( Solver % Values, & 1359 'Linear System Number of Threads', stat ) 1360 IF ( .NOT. stat ) nprocs = 1 1361! 1362 Factorize = ListGetLogical( Solver % Values, & 1363 'Linear System Refactorize', stat ) 1364 IF ( .NOT. stat ) Factorize = .TRUE. 1365 1366 IF ( Factorize .OR. A % SuperLU_Factors==0 ) THEN 1367 1368 IF ( A % SuperLU_Factors/= 0 ) THEN 1369 iopt = 3 1370 call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, & 1371 A % Rows, x, n, A % SuperLU_Factors, iinfo ) 1372 A % SuperLU_Factors=0 1373 END IF 1374 1375 ! First, factorize the matrix. The factors are stored in *factors* handle. 1376 iopt = 1 1377 call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, & 1378 A % Rows, x, n, A % SuperLU_Factors, iinfo ) 1379 1380 if (iinfo .eq. 0) then 1381 write (*,*) 'Factorization succeeded' 1382 else 1383 write(*,*) 'INFO from factorization = ', iinfo 1384 endif 1385 END IF 1386 1387 ! 1388 ! Second, solve the system using the existing factors. 1389 iopt = 2 1390 call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, & 1391 A % Rows, x, n, A % SuperLU_Factors, iinfo ) 1392! 1393 if (iinfo .eq. 0) then 1394 write (*,*) 'Solve succeeded' 1395 else 1396 write(*,*) 'INFO from triangular solve = ', iinfo 1397 endif 1398 1399! Last, free the storage allocated inside SuperLU 1400 FreeFactorize = ListGetLogical( Solver % Values, & 1401 'Linear System Free Factorization', stat ) 1402 IF ( .NOT. stat ) FreeFactorize = .TRUE. 1403 1404 IF ( Factorize .AND. FreeFactorize ) THEN 1405 iopt = 3 1406 call Solve_SuperLU( iopt, nprocs, n, nnz, nrhs, A % Values, A % Cols, & 1407 A % Rows, x, n, A % SuperLU_Factors, iinfo ) 1408 A % SuperLU_Factors = 0 1409 END IF 1410#endif 1411!------------------------------------------------------------------------------ 1412 END SUBROUTINE SuperLU_SolveSystem 1413!------------------------------------------------------------------------------ 1414 1415 1416!------------------------------------------------------------------------------ 1417!> Permon solver 1418!------------------------------------------------------------------------------ 1419 SUBROUTINE Permon_SolveSystem( Solver,A,x,b,Free_Fact ) 1420!------------------------------------------------------------------------------ 1421#ifdef HAVE_FETI4I 1422 use feti4i 1423#endif 1424 1425 LOGICAL, OPTIONAL :: Free_Fact 1426 TYPE(Matrix_t) :: A 1427 TYPE(Solver_t) :: Solver 1428 REAL(KIND=dp), TARGET :: x(*), b(*) 1429 1430#ifdef HAVE_FETI4I 1431 INCLUDE 'mpif.h' 1432 1433 INTEGER, ALLOCATABLE :: Owner(:) 1434 INTEGER :: i,j,n,nd,ip,ierr,icntlft,nzloc 1435 LOGICAL :: Factorize, FreeFactorize, stat, matsym, matspd, scaled 1436 1437 INTEGER :: n_dof_partition 1438 1439 INTEGER, ALLOCATABLE :: memb(:), DirichletInds(:), Neighbours(:), IntOptions(:) 1440 REAL(KIND=dp), ALLOCATABLE :: DirichletVals(:), RealOptions(:) 1441 INTEGER :: Comm_active, Group_active, Group_world 1442 1443 REAL(KIND=dp), ALLOCATABLE :: dbuf(:) 1444 1445 1446 n_dof_partition = A % NumberOfRows 1447 1448! INTERFACE 1449! FUNCTION Permon_InitSolve(n, gnum, nd, dinds, dvals, n_n, n_ranks) RESULT(handle) BIND(c,name='permon_initsolve') 1450! USE, INTRINSIC :: ISO_C_BINDING 1451! TYPE(C_PTR) :: handle 1452! INTEGER(C_INT), VALUE :: n, nd, n_n 1453! REAL(C_DOUBLE) :: dvals(*) 1454! INTEGER(C_INT) :: gnum(*), dinds(*), n_ranks(*) 1455! END FUNCTION Permon_Initsolve 1456! 1457! SUBROUTINE Permon_Solve( handle, x, b ) BIND(c,name='permon_solve') 1458! USE, INTRINSIC :: ISO_C_BINDING 1459! REAL(C_DOUBLE) :: x(*), b(*) 1460! TYPE(C_PTR), VALUE :: handle 1461! END SUBROUTINE Permon_solve 1462! END INTERFACE 1463 1464 IF ( PRESENT(Free_Fact) ) THEN 1465 IF ( Free_Fact ) THEN 1466 RETURN 1467 END IF 1468 END IF 1469 1470 Factorize = ListGetLogical( Solver % Values, 'Linear System Refactorize', stat ) 1471 IF ( .NOT. stat ) Factorize = .TRUE. 1472 1473 IF ( Factorize .OR. .NOT.C_ASSOCIATED(A % PermonSolverInstance) ) THEN 1474 IF ( C_ASSOCIATED(A % PermonSolverInstance) ) THEN 1475 CALL Fatal( 'Permon', 're-entry not implemented' ) 1476 END IF 1477 1478 nd = COUNT(A % ConstrainedDOF) 1479 ALLOCATE(DirichletInds(nd), DirichletVals(nd)) 1480 j = 0 1481 DO i=1,A % NumberOfRows 1482 IF(A % ConstrainedDOF(i)) THEN 1483 j = j + 1 1484 DirichletInds(j) = i; DirichletVals(j) = A % Dvalues(i) 1485 END IF 1486 END DO 1487 1488 !TODO sequential case not working 1489 n = 0 1490 ALLOCATE(neighbours(Parenv % PEs)) 1491 DO i=1,ParEnv % PEs 1492 IF( ParEnv % IsNeighbour(i) .AND. i-1/=ParEnv % myPE) THEN 1493 n = n + 1 1494 neighbours(n) = i-1 1495 END IF 1496 END DO 1497 1498 !A % PermonSolverInstance = Permon_InitSolve( SIZE(A % ParallelInfo % GlobalDOFs), & 1499 ! A % ParallelInfo % GlobalDOFs, nd, DirichletInds, DirichletVals, n, neighbours ) 1500 1501 IF( n_dof_partition /= SIZE(A % ParallelInfo % GlobalDOFs) ) THEN 1502 CALL Fatal( 'Permon', & 1503 'inconsistency: A % NumberOfRows /= SIZE(A % ParallelInfo % GlobalDOFs' ) 1504 END IF 1505 1506 ALLOCATE(IntOptions(10)) 1507 ALLOCATE(RealOptions(1)) 1508 1509 CALL FETI4ISetDefaultIntegerOptions(IntOptions) 1510 CALL FETI4ISetDefaultRealOptions(RealOptions) 1511 1512 CALL FETI4ICreateInstance(A % PermonSolverInstance, A % PermonMatrix, & 1513 A % NumberOfRows, b, A % ParallelInfo % GlobalDOFs, & 1514 n, neighbours, & 1515 nd, DirichletInds, DirichletVals, & 1516 IntOptions, RealOptions) 1517 END IF 1518 1519 !CALL Permon_Solve( A % PermonSolverInstance, x, b ) 1520 CALL FETI4ISolve(A % PermonSolverInstance, n_dof_partition, x) 1521#else 1522 CALL Fatal( 'Permon_SolveSystem', 'Permon Solver has not been installed.' ) 1523#endif 1524!------------------------------------------------------------------------------ 1525 END SUBROUTINE Permon_SolveSystem 1526!------------------------------------------------------------------------------ 1527 1528 1529 1530!------------------------------------------------------------------------------ 1531!> Solves a linear system using Pardiso direct solver (from either MKL or 1532!> official Pardiso distribution. If possible, MKL-version is used). 1533!------------------------------------------------------------------------------ 1534 SUBROUTINE Pardiso_SolveSystem( Solver,A,x,b,Free_fact ) 1535!------------------------------------------------------------------------------ 1536 IMPLICIT NONE 1537 1538 TYPE(Solver_t) :: Solver 1539 TYPE(Matrix_t) :: A 1540 REAL(KIND=dp), TARGET :: x(*), b(*) 1541 LOGICAL, OPTIONAL :: Free_fact 1542 1543! MKL version of Pardiso (interface is different) 1544#if defined(HAVE_MKL) 1545 INTERFACE 1546 SUBROUTINE pardiso(pt, maxfct, mnum, mtype, phase, n, & 1547 values, rows, cols, perm, nrhs, iparm, msglvl, b, x, ierror) 1548 USE Types 1549 IMPLICIT NONE 1550 REAL(KIND=dp) :: values(*), b(*), x(*) 1551 INTEGER(KIND=AddrInt) :: pt(*) 1552 INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror 1553 INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*) 1554 END SUBROUTINE pardiso 1555 1556 SUBROUTINE pardisoinit(pt, mtype, iparm) 1557 USE Types 1558 IMPLICIT NONE 1559 INTEGER(KIND=AddrInt) :: pt(*) 1560 INTEGER :: mtype 1561 INTEGER :: iparm(*) 1562 END SUBROUTINE pardisoinit 1563 END INTERFACE 1564 1565 INTEGER maxfct, mnum, mtype, phase, n, nrhs, ierror, msglvl 1566 INTEGER, POINTER :: Iparm(:) 1567 INTEGER i, j, k, nz, idum(1), nzutd 1568 LOGICAL :: Found, matsym, matpd 1569 REAL*8 :: ddum(1) 1570 1571 LOGICAL :: Factorize, FreeFactorize 1572 INTEGER :: tlen, allocstat 1573 CHARACTER(LEN=MAX_NAME_LEN) :: threads, mat_type 1574 1575 REAL(KIND=dp), POINTER :: values(:) 1576 INTEGER, POINTER :: rows(:), cols(:) 1577 1578 ! Check if system needs to be refactorized 1579 Factorize = ListGetLogical( Solver % Values, & 1580 'Linear System Refactorize', Found ) 1581 IF ( .NOT. Found ) Factorize = .TRUE. 1582 1583 ! Set matrix type for Pardiso 1584 mat_type = ListGetString( Solver % Values, 'Linear System Matrix Type', Found ) 1585 1586 IF (Found) THEN 1587 SELECT CASE(mat_type) 1588 CASE('positive definite') 1589 mtype = 2 1590 CASE('symmetric indefinite') 1591 mtype = -2 1592 CASE('structurally symmetric') 1593 mtype = 1 1594 CASE('nonsymmetric', 'general') 1595 mtype = 11 1596 CASE DEFAULT 1597 mtype = 11 1598 END SELECT 1599 ELSE 1600 ! Check if matrix is symmetric or spd 1601 matsym = ListGetLogical(Solver % Values, & 1602 'Linear System Symmetric', Found) 1603 1604 matpd = ListGetLogical(Solver % Values, & 1605 'Linear System Positive Definite', Found) 1606 1607 IF (matsym) THEN 1608 IF (matpd) THEN 1609 ! Matrix is symmetric positive definite 1610 mtype = 2 1611 ELSE 1612 ! Matrix is structurally symmetric (can't handle indefinite systems!!!!) 1613 mtype = 1 1614 END IF 1615 ELSE 1616 ! Matrix is unsymmetric 1617 mtype = 11 1618 END IF 1619 END IF 1620 1621 ! Free factorization if requested 1622 IF ( PRESENT(Free_Fact) ) THEN 1623 IF ( Free_Fact ) THEN 1624 IF(ASSOCIATED(A % PardisoId)) THEN 1625 phase = -1 ! release internal memory 1626 n = A % NumberOfRows 1627 maxfct = 1 1628 mnum = 1 1629 nrhs = 1 1630 msglvl = 0 1631 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1632 ddum, idum, idum, idum, nrhs, A % PardisoParam, msglvl, ddum, ddum, ierror) 1633 DEALLOCATE(A % PardisoId, A % PardisoParam) 1634 A % PardisoId => NULL() 1635 A % PardisoParam => NULL() 1636 END IF 1637 RETURN 1638 END IF 1639 END IF 1640 1641 ! Get number of rows and number of nonzero elements 1642 n = A % Numberofrows 1643 nz = A % Rows(n+1)-1 1644 1645 ! Copy upper triangular part of symmetric positive definite system 1646 IF ( ABS(mtype) == 2 ) THEN 1647 nzutd = 0 1648 DO i=1,n 1649 nzutd = nzutd + A % Rows(i+1)-A % Diag(i) 1650 END DO 1651 1652 ALLOCATE( values(nzutd), cols(nzutd), rows(n+1), STAT=allocstat) 1653 IF (allocstat /= 0) THEN 1654 CALL Fatal('Pardiso_SolveSystem', & 1655 'Memory allocation for row and column indices failed') 1656 END IF 1657 1658 ! Copy upper triangular part of A to Pardiso structure 1659 Rows(1)=1 1660 DO i=1,n 1661 ! Set up row pointers and copy values 1662 nzutd = A % Rows(i+1)-A % Diag(i) 1663 Rows(i+1) = Rows(i)+nzutd 1664 DO j=0,nzutd-1 1665 Cols(Rows(i)+j)=A % Cols(A%Diag(i)+j) 1666 Values(Rows(i)+j)=A % Values(A%Diag(i)+j) 1667 END DO 1668 END DO 1669 1670 ELSE 1671 Cols => A % Cols 1672 Rows => A % Rows 1673 Values => A % Values 1674 END IF 1675 1676 ! Set up parameters 1677 msglvl = 0 ! Do not write out any info 1678 maxfct = 1 ! Set up space for 1 matrix at most 1679 mnum = 1 ! Matrix to use in the solution phase (1st and only one) 1680 nrhs = 1 ! Use only one RHS 1681 iparm => A % PardisoParam 1682 1683 ! Compute factorization if necessary 1684 IF ( Factorize .OR. .NOT.ASSOCIATED(A % PardisoID) ) THEN 1685 ! Free factorization 1686 IF (ASSOCIATED(A % PardisoId)) THEN 1687 phase = -1 ! release internal memory 1688 CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, & 1689 ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror) 1690 DEALLOCATE(A % PardisoId, A % PardisoParam) 1691 A % PardisoId => NULL() 1692 A % PardisoParam => NULL() 1693 END IF 1694 1695 ! Allocate pardiso internal structures 1696 ALLOCATE(A % PardisoId(64), A % PardisoParam(64), STAT=allocstat) 1697 IF (allocstat /= 0) THEN 1698 CALL Fatal('Pardiso_SolveSystem', & 1699 'Memory allocation for Pardiso failed') 1700 END IF 1701 iparm => A % PardisoParam 1702 iparm = 0 1703 A % PardisoId = 0 1704 1705 ! Set up scaling values for solver based on matrix type 1706 CALL pardisoinit(A % PardisoId, mtype, iparm) 1707 1708 ! Set up rest of parameters explicitly 1709 iparm(1)=1 ! Do not use solver default parameters 1710 iparm(2)=2 ! Minimize fill-in with nested dissection from Metis 1711 iparm(3)=0 ! Reserved 1712 iparm(4)=0 ! Always compute factorization 1713 iparm(5)=0 ! No user input permutation 1714 iparm(6)=0 ! Write solution vector to x 1715 iparm(8)=5 ! Number of iterative refinement steps 1716 iparm(18)=-1 ! Report nnz(L) and nnz(U) 1717 iparm(19)=0 ! Do not report Mflops 1718 iparm(27)=0 ! Do not check sparse matrix representation 1719 iparm(35)=0 ! Use Fortran style indexing 1720 iparm(60)=0 ! Use in-core version of Pardiso 1721 1722 ! Perform analysis 1723 phase = 11 ! Analysis 1724 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1725 values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror) 1726 1727 IF (ierror .NE. 0) THEN 1728 WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror 1729 CALL Fatal('Pardiso_SolveSystem','Error during analysis phase') 1730 END IF 1731 1732 ! Perform factorization 1733 phase = 22 ! Factorization 1734 CALL pardiso (A % pardisoId, maxfct, mnum, mtype, phase, n, & 1735 values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror) 1736 1737 IF (ierror .NE. 0) THEN 1738 WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror 1739 CALL Fatal('Pardiso_SolveSystem','Error during factorization phase') 1740 END IF 1741 END IF ! Compute factorization 1742 1743 ! Perform solve 1744 phase = 33 ! Solve, iterative refinement 1745 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1746 values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror) 1747 1748 IF (ierror .NE. 0) THEN 1749 WRITE(*,'(A,I0)') 'MKL Pardiso: ERROR=', ierror 1750 CALL Fatal('Pardiso_SolveSystem','Error during solve phase') 1751 END IF 1752 1753 ! Release memory if needed 1754 FreeFactorize = ListGetLogical( Solver % Values, & 1755 'Linear System Free Factorization', Found ) 1756 IF ( .NOT. Found ) FreeFactorize = .TRUE. 1757 1758 IF ( Factorize .AND. FreeFactorize ) THEN 1759 phase = -1 ! release internal memory 1760 CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, & 1761 ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror) 1762 1763 DEALLOCATE(A % PardisoId, A % PardisoParam) 1764 A % PardisoId => NULL() 1765 A % PardisoParam => NULL() 1766 END IF 1767 IF (ABS(mtype) == 2) DEALLOCATE(Values, Rows, Cols) 1768 1769! Distribution version of Pardiso 1770#elif defined(HAVE_PARDISO) 1771!.. All other variables 1772 1773 INTERFACE 1774 SUBROUTINE pardiso(pt, maxfct, mnum, mtype, phase, n, & 1775 values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror, dparm) 1776 USE Types 1777 REAL(KIND=dp) :: values(*), b(*), x(*), dparm(*) 1778 INTEGER(KIND=AddrInt) :: pt(*) 1779 INTEGER :: idum(*), nrhs, iparm(*), msglvl, ierror 1780 INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*) 1781 END SUBROUTINE pardiso 1782 1783 SUBROUTINE pardisoinit(pt,mtype,solver,iparm,dparm,ierror) 1784 USE Types 1785 INTEGER :: mtype, iparm(*),ierror,solver 1786 REAL(KIND=dp) :: dparm(*) 1787 INTEGER(KIND=AddrInt) :: pt(*) 1788 END SUBROUTINE pardisoinit 1789 END INTERFACE 1790 1791 INTEGER maxfct, mnum, mtype, phase, n, nrhs, ierror, msglvl 1792 INTEGER, POINTER :: Iparm(:) 1793 INTEGER i, j, k, nz, idum(1) 1794 LOGICAL :: Found, Symm, Posdef 1795 REAL*8 waltime1, waltime2, ddum(1), dparm(64) 1796 1797 LOGICAL :: Factorize, FreeFactorize 1798 INTEGER :: tlen 1799 CHARACTER(LEN=MAX_STRING_LEN) :: threads 1800 1801 REAL(KIND=dp), POINTER :: values(:) 1802 INTEGER, POINTER :: rows(:), cols(:) 1803 1804 IF ( PRESENT(Free_Fact) ) THEN 1805 IF ( Free_Fact ) THEN 1806 IF(ASSOCIATED(A % PardisoId)) THEN 1807 phase = -1 ! release internal memory 1808 n = A % NumberOfRows 1809 mtype = 11 1810 maxfct = 1 1811 mnum = 1 1812 nrhs = 1 1813 msglvl = 0 1814 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1815 ddum, idum, idum, idum, nrhs, A % PardisoParam, msglvl, ddum, ddum, ierror,dparm) 1816 DEALLOCATE(A % PardisoId, A % PardisoParam) 1817 A % PardisoId => NULL() 1818 A % PardisoParam => NULL() 1819 END IF 1820 RETURN 1821 END IF 1822 END IF 1823 1824! .. Setup Pardiso control parameters und initialize the solvers 1825! internal address pointers. This is only necessary for the FIRST 1826! call of the PARDISO solver. 1827! 1828 Factorize = ListGetLogical( Solver % Values, & 1829 'Linear System Refactorize', Found ) 1830 IF ( .NOT. Found ) Factorize = .TRUE. 1831 1832 symm = ListGetLogical( Solver % Values, & 1833 'Linear System Symmetric', Found ) 1834 1835 posdef = ListGetLogical( Solver % Values, & 1836 'Linear System Positive Definite', Found ) 1837 1838 mtype = 11 1839 1840 cols => a % cols 1841 rows => a % rows 1842 values => a % values 1843 n = A % Numberofrows 1844 1845 IF ( posdef ) THEN 1846 nz = A % rows(n+1)-1 1847 allocate( values(nz), cols(nz), rows(n+1) ) 1848 k = 1 1849 do i=1,n 1850 rows(i)=k 1851 do j=a % rows(i), a % rows(i+1)-1 1852 if ( a % cols(j)>=i .and. a % values(j) /= 0._dp ) then 1853 cols(k) = a % cols(j) 1854 values(k) = a % values(j) 1855 k = k + 1 1856 end if 1857 end do 1858 end do 1859 rows(n+1)=k 1860 mtype = 2 1861 ELSE IF ( symm ) THEN 1862 mtype = 1 1863 END IF 1864 1865 msglvl = 0 ! with statistical information 1866 maxfct = 1 1867 mnum = 1 1868 nrhs = 1 1869 iparm => A % PardisoParam 1870 1871 1872 IF ( Factorize .OR. .NOT.ASSOCIATED(A % PardisoID) ) THEN 1873 IF(ASSOCIATED(A % PardisoId)) THEN 1874 phase = -1 ! release internal memory 1875 CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, & 1876 ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror,dparm) 1877 DEALLOCATE(A % PardisoId, A % PardisoParam) 1878 A % PardisoId => NULL() 1879 A % PardisoParam => NULL() 1880 END IF 1881 1882 ALLOCATE(A % PardisoId(64), A % PardisoParam(64)) 1883 iparm => A % PardisoParam 1884 CALL pardisoinit(A % PardisoId, mtype, 0, iparm, dparm, ierror ) 1885 1886!.. Reordering and Symbolic Factorization, This step also allocates 1887! all memory that is necessary for the factorization 1888! 1889 phase = 11 ! only reordering and symbolic factorization 1890 msglvl = 0 ! with statistical information 1891 maxfct = 1 1892 mnum = 1 1893 nrhs = 1 1894 1895! .. Numbers of Processors ( value of OMP_NUM_THREADS ) 1896 CALL envir( 'OMP_NUM_THREADS', threads, tlen ) 1897 1898 iparm(3) = 1 1899 IF ( tlen>0 ) & 1900 READ(threads(1:tlen),*) iparm(3) 1901 IF (iparm(3)<=0) Iparm(3) = 1 1902 1903 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1904 values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm) 1905 1906 IF (ierror .NE. 0) THEN 1907 WRITE(*,*) 'The following ERROR was detected: ', ierror 1908 STOP EXIT_ERROR 1909 END IF 1910 1911!.. Factorization. 1912 phase = 22 ! only factorization 1913 CALL pardiso (A % pardisoId, maxfct, mnum, mtype, phase, n, & 1914 values, rows, cols, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm) 1915 1916 IF (ierror .NE. 0) THEN 1917 WRITE(*,*) 'The following ERROR was detected: ', ierror 1918 STOP EXIT_ERROR 1919 ENDIF 1920 END IF 1921 1922!.. Back substitution and iterative refinement 1923 phase = 33 ! only factorization 1924 iparm(8) = 0 ! max number of iterative refinement steps 1925 ! (if perturbations occur in the factorization 1926 ! phase, two refinement steps are taken) 1927 1928 CALL pardiso(A % PardisoId, maxfct, mnum, mtype, phase, n, & 1929 values, rows, cols, idum, nrhs, iparm, msglvl, b, x, ierror, dparm) 1930 1931!.. Termination and release of memory 1932 FreeFactorize = ListGetLogical( Solver % Values, & 1933 'Linear System Free Factorization', Found ) 1934 IF ( .NOT. Found ) FreeFactorize = .TRUE. 1935 1936 IF ( Factorize .AND. FreeFactorize ) THEN 1937 phase = -1 ! release internal memory 1938 CALL pardiso (A % PardisoId, maxfct, mnum, mtype, phase, n, & 1939 ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, ierror, dparm) 1940 1941 DEALLOCATE(A % PardisoId, A % PardisoParam) 1942 A % PardisoId => NULL() 1943 A % PardisoParam => NULL() 1944 END IF 1945 1946 IF(posdef) DEALLOCATE( values, rows, cols ) 1947#else 1948 CALL Fatal( 'Parsido_SolveSystem', 'Pardiso solver has not been installed.' ) 1949#endif 1950 1951!------------------------------------------------------------------------------ 1952 END SUBROUTINE Pardiso_SolveSystem 1953!------------------------------------------------------------------------------ 1954 1955!------------------------------------------------------------------------------ 1956!> Solves a linear system using Cluster Pardiso direct solver from MKL 1957!------------------------------------------------------------------------------ 1958 SUBROUTINE CPardiso_SolveSystem( Solver,A,x,b,Free_fact ) 1959!------------------------------------------------------------------------------ 1960 IMPLICIT NONE 1961 1962 TYPE(Solver_t) :: Solver 1963 TYPE(Matrix_t) :: A 1964 REAL(KIND=dp), TARGET :: x(*), b(*) 1965 LOGICAL, OPTIONAL :: Free_fact 1966 1967! Cluster Pardiso 1968#if defined(HAVE_MKL) && defined(HAVE_CPARDISO) 1969 INTERFACE 1970 SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, & 1971 values, rows, cols, perm, nrhs, iparm, msglvl, b, x, comm, ierror) 1972 USE Types 1973 REAL(KIND=dp) :: values(*), b(*), x(*) 1974 INTEGER(KIND=AddrInt) :: pt(*) 1975 INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror 1976 INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm 1977 END SUBROUTINE cluster_sparse_solver 1978 END INTERFACE 1979 1980 INTEGER :: phase, n, ierror 1981 INTEGER, POINTER :: Iparm(:) 1982 INTEGER i, j, k, nz, nzutd, idum(1), nl, nt 1983 LOGICAL :: Found, matsym, matpd 1984 REAL(kind=dp) :: ddum(1) 1985 REAL(kind=dp), POINTER, DIMENSION(:) :: dbuf 1986 1987 LOGICAL :: Factorize, FreeFactorize 1988 INTEGER :: tlen, allocstat 1989 1990 REAL(KIND=dp), POINTER CONTIG :: values(:) 1991 INTEGER, POINTER CONTIG :: rows(:), cols(:) 1992 1993 ! Free factorization if requested 1994 IF ( PRESENT(Free_Fact) ) THEN 1995 IF ( Free_Fact ) THEN 1996 CALL CPardiso_Free(A) 1997 END IF 1998 1999 RETURN 2000 END IF 2001 2002 ! Check if system needs to be refactorized 2003 Factorize = ListGetLogical( Solver % Values, & 2004 'Linear System Refactorize', Found ) 2005 IF ( .NOT. Found ) Factorize = .TRUE. 2006 2007 ! Compute factorization if necessary 2008 IF ( Factorize .OR. .NOT.ASSOCIATED(A % CPardisoID) ) THEN 2009 CALL CPardiso_Factorize(Solver, A) 2010 END IF 2011 2012 ! Get global start and end of domain 2013 nl = A % CPardisoId % iparm(41) 2014 nt = A % CPardisoId % iparm(42) 2015 2016 ! Gather RHS 2017 A % CPardisoId % rhs = 0D0 2018 DO i=1,A % NumberOfRows 2019 A % CPardisoId % rhs(A % Gorder(i)-nl+1) = b(i) 2020 END DO 2021 2022 ! Perform solve 2023 phase = 33 ! Solve, iterative refinement 2024 CALL cluster_sparse_solver(A % CPardisoId % ID, & 2025 A % CPardisoId % maxfct, A % CPardisoId % mnum, & 2026 A % CPardisoId % mtype, phase, A % CPardisoId % n, & 2027 A % CPardisoID % aa, A % CPardisoID % ia, A % CPardisoID % ja, idum, & 2028 A % CPardisoId % nrhs, A % CPardisoID % iparm, & 2029 A % CPardisoId % msglvl, A % CPardisoId % rhs, & 2030 A % CPardisoId % x, A % Comm, ierror) 2031 2032 IF (ierror /= 0) THEN 2033 WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror 2034 CALL Fatal('CPardiso_SolveSystem','Error during solve phase') 2035 END IF 2036 2037 ! Distribute solution 2038 DO i=1,A % NumberOfRows 2039 x(i)=A % CPardisoId % x(A % Gorder(i)-nl+1) 2040 END DO 2041 2042 ! Release memory if needed 2043 FreeFactorize = ListGetLogical( Solver % Values, & 2044 'Linear System Free Factorization', Found ) 2045 IF ( .NOT. Found ) FreeFactorize = .TRUE. 2046 2047 IF ( Factorize .AND. FreeFactorize ) THEN 2048 CALL CPardiso_Free(A) 2049 END IF 2050 2051#else 2052 CALL Fatal( 'CParsido_SolveSystem', 'Cluster Pardiso solver has not been installed.' ) 2053#endif 2054!------------------------------------------------------------------------------ 2055 END SUBROUTINE CPardiso_SolveSystem 2056!------------------------------------------------------------------------------ 2057 2058#if defined(HAVE_MKL) && defined(HAVE_CPARDISO) 2059 SUBROUTINE CPardiso_Factorize(Solver, A) 2060 IMPLICIT NONE 2061 TYPE(Solver_t) :: Solver 2062 TYPE(Matrix_t) :: A 2063 2064 INTERFACE 2065 SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, & 2066 values, rows, cols, perm, nrhs, iparm, msglvl, b, x, comm, ierror) 2067 USE Types 2068 REAL(KIND=dp) :: values(*), b(*), x(*) 2069 INTEGER(KIND=AddrInt) :: pt(*) 2070 INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror 2071 INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm 2072 END SUBROUTINE cluster_sparse_solver 2073 END INTERFACE 2074 2075 LOGICAL :: matsym, matpd, Found 2076 INTEGER :: i, j, k, rind, lrow, rptr, rsize, lind, tind 2077 INTEGER :: allocstat, n, nz, nzutd, nl, nt, nOwned, nhalo, ierror 2078 INTEGER :: phase, idum(1) 2079 REAL(kind=dp) :: ddum(1) 2080 REAL(kind=dp), DIMENSION(:), POINTER :: aa 2081 INTEGER, DIMENSION(:), POINTER CONTIG :: iparm, ia, ja, owner, dsize, iperm, Order 2082 2083 INTEGER :: fid 2084 CHARACTER(LEN=MAX_NAME_LEN) :: mat_type 2085 2086 ! Free old factorization if necessary 2087 IF (ASSOCIATED(A % CPardisoId)) THEN 2088 CALL CPardiso_Free(A) 2089 END IF 2090 2091 ! Allocate Pardiso structure 2092 ALLOCATE(A % CPardisoId, STAT=allocstat) 2093 IF (allocstat /= 0) THEN 2094 CALL Fatal('CPardiso_Factorize', & 2095 'Memory allocation for CPardiso failed') 2096 END IF 2097 ! Allocate control structures 2098 ALLOCATE(A % CPardisoId% ID(64), A % CPardisoId % IParm(64), STAT=allocstat) 2099 IF (allocstat /= 0) THEN 2100 CALL Fatal('CPardiso_Factorize', & 2101 'Memory allocation for CPardiso failed') 2102 END IF 2103 2104 iparm => A % CPardisoId % IParm 2105 ! Initialize control parameters and solver Id's 2106 DO i=1,64 2107 iparm(i)=0 2108 END DO 2109 DO i=1,64 2110 A % CPardisoId % ID(i)=0 2111 END DO 2112 2113 ! Set matrix type for CPardiso 2114 mat_type = ListGetString( Solver % Values, 'Linear System Matrix Type', Found ) 2115 IF (Found) THEN 2116 SELECT CASE(mat_type) 2117 CASE('positive definite') 2118 A % CPardisoID % mtype = 2 2119 CASE('symmetric indefinite') 2120 A % CPardisoID % mtype = -2 2121 CASE('structurally symmetric') 2122 A % CPardisoID % mtype = 1 2123 CASE('nonsymmetric', 'general') 2124 A % CPardisoID % mtype = 11 2125 CASE DEFAULT 2126 A % CPardisoID % mtype = 11 2127 END SELECT 2128 ELSE 2129 ! Check if matrix is symmetric or spd 2130 matsym = ListGetLogical(Solver % Values, & 2131 'Linear System Symmetric', Found) 2132 2133 matpd = ListGetLogical(Solver % Values, & 2134 'Linear System Positive Definite', Found) 2135 2136 IF (matsym) THEN 2137 IF (matpd) THEN 2138 ! Matrix is symmetric positive definite 2139 A % CPardisoID % mtype = 2 2140 ELSE 2141 ! Matrix is structurally symmetric 2142 A % CPardisoID % mtype = 1 2143 END IF 2144 ELSE 2145 ! Matrix is nonsymmetric 2146 A % CPardisoID % mtype = 11 2147 END IF 2148 END IF 2149 2150 ! Set up continuous numbering for the whole computation domain 2151 n = SIZE(A % ParallelInfo % GlobalDOFs) 2152 ALLOCATE(A % Gorder(n), Owner(n), STAT=allocstat) 2153 IF (allocstat /= 0) THEN 2154 CALL Fatal('CPardiso_Factorize', & 2155 'Memory allocation for CPardiso global numbering failed') 2156 END IF 2157 CALL ContinuousNumbering(A % ParallelInfo, A % Perm, A % Gorder, Owner, nOwn=nOwned) 2158 2159 ! Compute the number of global dofs 2160 CALL MPI_ALLREDUCE(nOwned, A % CPardisoId % n, & 2161 1, MPI_INTEGER, MPI_SUM, A % Comm, ierror) 2162 DEALLOCATE(Owner) 2163 2164 ! Find bounds of domain 2165 nl = A % Gorder(1) 2166 nt = A % Gorder(1) 2167 DO i=2,n 2168 ! NOTE: Matrix is structurally symmetric 2169 rind = A % Gorder(i) 2170 nl = MIN(rind, nl) 2171 nt = MAX(rind, nt) 2172 END DO 2173 2174 ! Allocate temp storage for global numbering 2175 ALLOCATE(Order(n), iperm(n), STAT=allocstat) 2176 IF (allocstat /= 0) THEN 2177 CALL Fatal('CPardiso_Factorize', & 2178 'Memory allocation for CPardiso global numbering failed') 2179 END IF 2180 2181 ! Sort global numbering to build matrix 2182 Order(1:n) = A % Gorder(1:n) 2183 DO i=1,n 2184 iperm(i)=i 2185 END DO 2186 CALL SortI(n, Order, iperm) 2187 2188 ! Allocate storage for CPardiso matrix 2189 nhalo = (nt-nl+1)-n 2190 nz = A % Rows(A % NumberOfRows+1)-1 2191 IF (ABS(A % CPardisoID % mtype) == 2) THEN 2192 nzutd = ((nz-n)/2)+1 + n 2193 ALLOCATE(A % CPardisoID % ia(nt-nl+2), & 2194 A % CPardisoId % ja(nzutd+nhalo), & 2195 A % CPardisoId % aa(nzutd+nhalo), & 2196 STAT=allocstat) 2197 ELSE 2198 ALLOCATE(A % CPardisoID % ia(nt-nl+2), & 2199 A % CPardisoId % ja(nz+nhalo), & 2200 A % CPardisoId % aa(nz+nhalo), & 2201 STAT=allocstat) 2202 END IF 2203 IF (allocstat /= 0) THEN 2204 CALL Fatal('CPardiso_Factorize', & 2205 'Memory allocation for CPardiso matrix failed') 2206 END IF 2207 ia => A % CPardisoId % ia 2208 ja => A % CPardisoId % ja 2209 aa => A % CPardisoId % aa 2210 2211 ! Build distributed CRS matrix 2212 ia(1) = 1 2213 lrow = 1 ! Next row to add 2214 rptr = 1 ! Pointer to next row to add, equals ia(lrow) 2215 lind = Order(1)-1 ! Row pointer for the first round 2216 2217 ! Add rows of matrix 2218 DO i=1,n 2219 ! Add empty rows until the beginning of the row to add 2220 ! (first round adds nothing due to choice of lind) 2221 tind = Order(i) 2222 rsize = (tind-lind)-1 2223 2224 ! Put zeroes to the diagonal 2225 DO j=1,rsize 2226 ia(lrow+j)=rptr+j 2227 ja(rptr+(j-1))=lind+j 2228 aa(rptr+(j-1))=0D0 2229 END DO 2230 ! Set up row pointers 2231 rptr = rptr + rsize 2232 lrow = lrow + rsize 2233 2234 ! Add next row 2235 rind = iperm(i) 2236 lind = A % rows(rind) 2237 tind = A % rows(rind+1) 2238 IF (ABS(A % CPardisoId % mtype) == 2) THEN 2239 ! Add only upper triangular elements for symmetric matrices 2240 rsize = 0 2241 DO j=lind, tind-1 2242 IF (A % Gorder(A % Cols(j)) >= Order(i)) THEN 2243 ja(rptr+rsize)=A % Gorder(A % Cols(j)) 2244 aa(rptr+rsize)=A % values(j) 2245 rsize = rsize + 1 2246 END IF 2247 END DO 2248 2249 ELSE 2250 rsize = tind-lind 2251 DO j=lind, tind-1 2252 ja(rptr+(j-lind))=A % Gorder(A % Cols(j)) 2253 aa(rptr+(j-lind))=A % values(j) 2254 END DO 2255 END IF 2256 2257 ! Sort column indices 2258 CALL SortF(rsize, ja(rptr:rptr+rsize), aa(rptr:rptr+rsize)) 2259 2260 ! Set up row pointers 2261 rptr = rptr + rsize 2262 lrow = lrow + 1 2263 ia(lrow) = rptr 2264 2265 lind = Order(i) ! Store row index for next round 2266 END DO 2267 2268 ! Deallocate temp storage 2269 DEALLOCATE(Order, iperm) 2270 2271 ! Set up parameters 2272 A % CPardisoId % msglvl = 0 ! Do not write out = 0 / write out = 1 info 2273 A % CPardisoId % maxfct = 1 ! Set up space for 1 matrix at most 2274 A % CPardisoId % mnum = 1 ! Matrix to use in the solution phase (1st and only one) 2275 A % CPardisoId % nrhs = 1 ! Use only one RHS 2276 ALLOCATE(A % CPardisoId % rhs(nt-nl+1), & 2277 A % CPardisoId % x(nt-nl+1), STAT=allocstat) 2278 IF (allocstat /= 0) THEN 2279 CALL Fatal('CPardiso_Factorize', & 2280 'Memory allocation for CPardiso rhs and solution vector x failed') 2281 END IF 2282 2283 ! Set up parameters explicitly 2284 iparm(1)=1 ! Do not use = 1 / use = 0 solver default parameters 2285 iparm(2)=2 ! Minimize fill-in with OpenMP nested dissection 2286 iparm(5)=0 ! No user input permutation 2287 iparm(6)=0 ! Write solution vector to x 2288 iparm(8)=0 ! Number of iterative refinement steps 2289 IF (A % CPardisoID % mtype ==11 .OR. & 2290 A % CPardisoID % mtype == 13) THEN 2291 iparm(10)=13 ! Perturbation value 10^-iparm(10) in case of small pivots 2292 iparm(11)=1 ! Use scalings from symmetric weighted matching 2293 iparm(13)=1 ! Use permutations from nonsymmetric weighted matching 2294 ELSE 2295 iparm(10)=8 ! Perturbation value 10^-iparm(10) in case of small pivots 2296 iparm(11)=0 ! Do not use scalings from symmetric weighted matching 2297 iparm(13)=0 ! Do not use permutations from symmetric weighted matching 2298 END IF 2299 2300 iparm(21)=1 ! Do not use Bunch Kaufman pivoting 2301 iparm(27)=0 ! Do not check sparse matrix representation 2302 iparm(28)=0 ! Use double precision 2303 iparm(35)=0 ! Use Fortran indexing 2304 2305 ! CPardiso matrix input format 2306 iparm(40) = 2 ! Distributed solution phase, distributed solution vector 2307 iparm(41) = nl ! Beginning of solution domain 2308 iparm(42) = nt ! End of solution domain 2309 2310 ! Perform analysis 2311 phase = 11 ! Analysis 2312 CALL cluster_sparse_solver(A % CPardisoId % ID, & 2313 A % CPardisoId % maxfct, A % CPardisoId % mnum, & 2314 A % CPardisoId % mtype, phase, A % CPardisoId % n, & 2315 aa, ia, ja, idum, A % CPardisoId % nrhs, iparm, & 2316 A % CPardisoId % msglvl, & 2317 ddum, ddum, A % Comm, ierror) 2318 IF (ierror /= 0) THEN 2319 WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror 2320 CALL Fatal('CPardiso_SolveSystem','Error during analysis phase') 2321 END IF 2322 2323 ! Perform factorization 2324 phase = 22 ! Factorization 2325 CALL cluster_sparse_solver(A % CPardisoId % ID, & 2326 A % CPardisoId % maxfct, A % CPardisoId % mnum, & 2327 A % CPardisoId % mtype, phase, A % CPardisoId % n, & 2328 aa, ia, ja, idum, A % CPardisoId % nrhs, iparm, & 2329 A % CPardisoId % msglvl, & 2330 ddum, ddum, A % Comm, ierror) 2331 IF (ierror .NE. 0) THEN 2332 WRITE(*,'(A,I0)') 'MKL CPardiso: ERROR=', ierror 2333 CALL Fatal('CPardiso_SolveSystem','Error during factorization phase') 2334 END IF 2335 END SUBROUTINE CPardiso_Factorize 2336 2337 2338 SUBROUTINE CPardiso_Free(A) 2339 IMPLICIT NONE 2340 2341 TYPE(Matrix_t) :: A 2342 INTERFACE 2343 SUBROUTINE cluster_sparse_solver(pt, maxfct, mnum, mtype, phase, n, & 2344 values, rows, cols, perm, nrhs, iparm, & 2345 msglvl, b, x, comm, ierror) 2346 USE Types 2347 REAL(KIND=dp) :: values(*), b(*), x(*) 2348 INTEGER(KIND=AddrInt) :: pt(*) 2349 INTEGER :: perm(*), nrhs, iparm(*), msglvl, ierror 2350 INTEGER :: maxfct, mnum, mtype, phase, n, rows(*), cols(*), comm 2351 END SUBROUTINE cluster_sparse_solver 2352 END INTERFACE 2353 2354 INTEGER :: ierror, phase, idum(1) 2355 REAL(kind=dp) :: ddum(1) 2356 2357 IF(ASSOCIATED(A % CPardisoId)) THEN 2358 phase = -1 ! release internal memory 2359 CALL cluster_sparse_solver(A % CPardisoId % ID, & 2360 A % CPardisoID % maxfct, A % CPardisoID % mnum, & 2361 A % CPardisoID % mtype, phase, A % CPardisoID % n, & 2362 A % CPardisoId % aa, A % CPardisoId % ia, A % CPardisoId % ja, & 2363 idum, A % CPardisoID % nrhs, A % CPardisoID % IParm, & 2364 A % CPardisoID % msglvl, ddum, ddum, A % Comm, ierror) 2365 2366 ! Deallocate global ordering 2367 DEALLOCATE(A % Gorder) 2368 2369 ! Deallocate control structures 2370 DEALLOCATE(A % CPardisoId % ID) 2371 DEALLOCATE(A % CPardisoID % IParm) 2372 ! Deallocate matrix and permutation 2373 DEALLOCATE(A % CPardisoId % ia, A % CPardisoId % ja, & 2374 A % CPardisoId % aa, A % CPardisoID % rhs, & 2375 A % CPardisoID % x) 2376 ! Deallocate CPardiso structure 2377 DEALLOCATE(A % CPardisoID) 2378 END IF 2379 END SUBROUTINE CPardiso_Free 2380#endif 2381 2382 2383!------------------------------------------------------------------------------ 2384 SUBROUTINE DirectSolver( A,x,b,Solver,Free_Fact ) 2385!------------------------------------------------------------------------------ 2386 2387 TYPE(Solver_t) :: Solver 2388 REAL(KIND=dp) :: x(*),b(*) 2389 TYPE(Matrix_t) :: A 2390 LOGICAL, OPTIONAL :: Free_Fact 2391!------------------------------------------------------------------------------ 2392 2393 LOGICAL :: GotIt 2394 CHARACTER(LEN=MAX_NAME_LEN) :: Method 2395!------------------------------------------------------------------------------ 2396 2397 IF ( PRESENT(Free_Fact) ) THEN 2398 IF ( Free_Fact ) THEN 2399 CALL BandSolver( A, x, b, Free_Fact ) 2400 CALL ComplexBandSolver( A, x, b, Free_Fact ) 2401#ifdef HAVE_MUMPS 2402 CALL Mumps_SolveSystem( Solver, A, x, b, Free_Fact ) 2403 CALL MumpsLocal_SolveSystem( Solver, A, x, b, Free_Fact ) 2404#endif 2405#if defined(HAVE_MKL) || defined(HAVE_PARDISO) 2406 CALL Pardiso_SolveSystem( Solver, A, x, b, Free_Fact ) 2407#endif 2408#if defined(HAVE_MKL) && defined(HAVE_CPARDISO) 2409 CALL CPardiso_SolveSystem( Solver, A, x, b, Free_Fact ) 2410#endif 2411#ifdef HAVE_SUPERLU 2412 CALL SuperLU_SolveSystem( Solver, A, x, b, Free_Fact ) 2413#endif 2414#ifdef HAVE_UMFPACK 2415 CALL Umfpack_SolveSystem( Solver, A, x, b, Free_Fact ) 2416#endif 2417#ifdef HAVE_CHOLMOD 2418 CALL SPQR_SolveSystem( Solver, A, x, b, Free_Fact ) 2419 CALL Cholmod_SolveSystem( Solver, A, x, b, Free_Fact ) 2420#endif 2421#ifdef HAVE_FETI4I 2422 CALL Permon_SolveSystem( Solver, A, x, b, Free_Fact ) 2423#endif 2424 RETURN 2425 RETURN 2426 END IF 2427 END IF 2428 2429 Method = ListGetString( Solver % Values, 'Linear System Direct Method',GotIt ) 2430 IF ( .NOT. GotIt ) Method = 'banded' 2431 2432 2433 CALL Info('DirectSolver','Using direct method: '//TRIM(Method),Level=9) 2434 2435 SELECT CASE( Method ) 2436 CASE( 'banded', 'symmetric banded' ) 2437 IF ( .NOT. A % Complex ) THEN 2438 CALL BandSolver( A, x, b ) 2439 ELSE 2440 CALL ComplexBandSolver( A, x, b ) 2441 END IF 2442 2443 CASE( 'umfpack', 'big umfpack' ) 2444 CALL Umfpack_SolveSystem( Solver, A, x, b ) 2445 2446 CASE( 'cholmod' ) 2447 CALL Cholmod_SolveSystem( Solver, A, x, b ) 2448 2449 CASE( 'spqr' ) 2450 CALL SPQR_SolveSystem( Solver, A, x, b ) 2451 2452 CASE( 'mumps' ) 2453 CALL Mumps_SolveSystem( Solver, A, x, b ) 2454 2455 CASE( 'mumpslocal' ) 2456 CALL MumpsLocal_SolveSystem( Solver, A, x, b ) 2457 2458 CASE( 'superlu' ) 2459 CALL SuperLU_SolveSystem( Solver, A, x, b ) 2460 2461 CASE( 'permon' ) 2462 CALL Permon_SolveSystem( Solver, A, x, b ) 2463 2464 CASE( 'pardiso' ) 2465 CALL Pardiso_SolveSystem( Solver, A, x, b ) 2466 2467 CASE( 'cpardiso' ) 2468 CALL CPardiso_SolveSystem( Solver, A, x, b ) 2469 2470 CASE DEFAULT 2471 CALL Fatal( 'DirectSolver', 'Unknown direct solver method.' ) 2472 END SELECT 2473 2474 ! We should be able to trust that a direct strategy will return a converged 2475 ! linear system. 2476 IF( ASSOCIATED( Solver % Variable ) ) THEN 2477 Solver % Variable % LinConverged = 1 2478 END IF 2479 2480!------------------------------------------------------------------------------ 2481 END SUBROUTINE DirectSolver 2482!------------------------------------------------------------------------------ 2483 2484END MODULE DirectSolve 2485 2486!> \} ElmerLib 2487