1! Copyright 1998-2019 Lawrence Livermore National Security, LLC and other 2! HYPRE Project Developers. See the top-level COPYRIGHT file for details. 3! 4! SPDX-License-Identifier: (Apache-2.0 OR MIT) 5 6! 7! Example 5 8! 9! Interface: Linear-Algebraic (IJ), Fortran (77) version 10! 11! Compile with: make ex5f 12! 13! Sample run: mpirun -np 4 ex5f 14! 15! Description: This example solves the 2-D 16! Laplacian problem with zero boundary conditions 17! on an nxn grid. The number of unknowns is N=n^2. 18! The standard 5-point stencil is used, and we solve 19! for the interior nodes only. 20! 21! This example solves the same problem as Example 3. 22! Available solvers are AMG, PCG, and PCG with AMG, 23! and PCG with ParaSails 24! 25! 26! Notes: for PCG, GMRES and BiCGStab, precond_id means: 27! 0 - do not set up a preconditioner 28! 1 - set up a ds preconditioner 29! 2 - set up an amg preconditioner 30! 3 - set up a pilut preconditioner 31! 4 - set up a ParaSails preconditioner 32! 33 34 program ex5f 35 36 37 implicit none 38 39 include 'mpif.h' 40 41 integer MAX_LOCAL_SIZE 42 integer HYPRE_PARCSR 43 44 parameter (MAX_LOCAL_SIZE=123000) 45 46! the following is from HYPRE.c 47 parameter (HYPRE_PARCSR=5555) 48 49 integer ierr 50 integer num_procs, myid 51 integer local_size, extra 52 integer n, solver_id, print_solution, ng 53 integer nnz, ilower, iupper, i 54 integer precond_id; 55 double precision h, h2 56 double precision rhs_values(MAX_LOCAL_SIZE) 57 double precision x_values(MAX_LOCAL_SIZE) 58 integer rows(MAX_LOCAL_SIZE) 59 integer cols(5) 60 double precision values(5) 61 integer num_iterations 62 double precision final_res_norm, tol 63 64 integer mpi_comm 65 66 integer*8 parcsr_A 67 integer*8 A 68 integer*8 b 69 integer*8 x 70 integer*8 par_b 71 integer*8 par_x 72 integer*8 solver 73 integer*8 precond 74 75!----------------------------------------------------------------------- 76! Initialize MPI 77!----------------------------------------------------------------------- 78 79 call MPI_INIT(ierr) 80 call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr) 81 call MPI_COMM_SIZE(MPI_COMM_WORLD, num_procs, ierr) 82 mpi_comm = MPI_COMM_WORLD 83 84 call HYPRE_Init(ierr) 85 86! Default problem parameters 87 n = 33 88 solver_id = 0 89 print_solution = 0 90 tol = 1.0d-7 91 92! The input section not implemented yet. 93 94! Preliminaries: want at least one processor per row 95 if ( n*n .lt. num_procs ) then 96 n = int(sqrt(real(num_procs))) + 1 97 endif 98! ng = global no. rows, h = mesh size 99 ng = n*n 100 h = 1.0d0/(n+1) 101 h2 = h*h 102 103! Each processor knows only of its own rows - the range is denoted by ilower 104! and upper. Here we partition the rows. We account for the fact that 105! N may not divide evenly by the number of processors. 106 local_size = ng/num_procs 107 extra = ng - local_size*num_procs 108 109 ilower = local_size*myid 110 ilower = ilower + min(myid, extra) 111 112 iupper = local_size*(myid+1) 113 iupper = iupper + min(myid+1, extra) 114 iupper = iupper - 1 115 116! How many rows do I have? 117 local_size = iupper - ilower + 1 118 119! Create the matrix. 120! Note that this is a square matrix, so we indicate the row partition 121! size twice (since number of rows = number of cols) 122 call HYPRE_IJMatrixCreate(mpi_comm, ilower, 123 1 iupper, ilower, iupper, A, ierr) 124 125 126! Choose a parallel csr format storage (see the User's Manual) 127 call HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR, ierr) 128 129! Initialize before setting coefficients 130 call HYPRE_IJMatrixInitialize(A, ierr) 131 132 133! Now go through my local rows and set the matrix entries. 134! Each row has at most 5 entries. For example, if n=3: 135! 136! A = [M -I 0; -I M -I; 0 -I M] 137! M = [4 -1 0; -1 4 -1; 0 -1 4] 138! 139! Note that here we are setting one row at a time, though 140! one could set all the rows together (see the User's Manual). 141 142 143 do i = ilower, iupper 144 nnz = 1 145 146 147! The left identity block:position i-n 148 if ( (i-n) .ge. 0 ) then 149 cols(nnz) = i-n 150 values(nnz) = -1.0d0 151 nnz = nnz + 1 152 endif 153 154! The left -1: position i-1 155 if ( mod(i,n).ne.0 ) then 156 cols(nnz) = i-1 157 values(nnz) = -1.0d0 158 nnz = nnz + 1 159 endif 160 161! Set the diagonal: position i 162 cols(nnz) = i 163 values(nnz) = 4.0d0 164 nnz = nnz + 1 165 166! The right -1: position i+1 167 if ( mod((i+1),n) .ne. 0 ) then 168 cols(nnz) = i+1 169 values(nnz) = -1.0d0 170 nnz = nnz + 1 171 endif 172 173! The right identity block:position i+n 174 if ( (i+n) .lt. ng ) then 175 cols(nnz) = i+n 176 values(nnz) = -1.0d0 177 nnz = nnz + 1 178 endif 179 180! Set the values for row i 181 call HYPRE_IJMatrixSetValues( 182 1 A, 1, nnz-1, i, cols, values, ierr) 183 184 enddo 185 186 187! Assemble after setting the coefficients 188 call HYPRE_IJMatrixAssemble(A, ierr) 189 190! Get parcsr matrix object 191 call HYPRE_IJMatrixGetObject(A, parcsr_A, ierr) 192 193 194! Create the rhs and solution 195 call HYPRE_IJVectorCreate(mpi_comm, 196 1 ilower, iupper, b, ierr) 197 call HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR, ierr) 198 call HYPRE_IJVectorInitialize(b, ierr) 199 200 call HYPRE_IJVectorCreate(mpi_comm, 201 1 ilower, iupper, x, ierr) 202 call HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR, ierr) 203 call HYPRE_IJVectorInitialize(x, ierr) 204 205 206! Set the rhs values to h^2 and the solution to zero 207 do i = 1, local_size 208 rhs_values(i) = h2 209 x_values(i) = 0.0 210 rows(i) = ilower + i -1 211 enddo 212 call HYPRE_IJVectorSetValues( 213 1 b, local_size, rows, rhs_values, ierr) 214 call HYPRE_IJVectorSetValues( 215 1 x, local_size, rows, x_values, ierr) 216 217 218 call HYPRE_IJVectorAssemble(b, ierr) 219 call HYPRE_IJVectorAssemble(x, ierr) 220 221! get the x and b objects 222 223 call HYPRE_IJVectorGetObject(b, par_b, ierr) 224 call HYPRE_IJVectorGetObject(x, par_x, ierr) 225 226 227! Choose a solver and solve the system 228 229! AMG 230 if ( solver_id .eq. 0 ) then 231 232! Create solver 233 call HYPRE_BoomerAMGCreate(solver, ierr) 234 235 236! Set some parameters (See Reference Manual for more parameters) 237 238! print solve info + parameters 239 call HYPRE_BoomerAMGSetPrintLevel(solver, 3, ierr) 240! old defaults, Falgout coarsening, mod. class. interpolation 241 call HYPRE_BoomerAMGSetOldDefault(solver, ierr) 242! G-S/Jacobi hybrid relaxation 243 call HYPRE_BoomerAMGSetRelaxType(solver, 3, ierr) 244! C/F relaxation 245 call HYPRE_BoomerAMGSetRelaxOrder(solver, 1, ierr) 246! Sweeeps on each level 247 call HYPRE_BoomerAMGSetNumSweeps(solver, 1, ierr) 248! maximum number of levels 249 call HYPRE_BoomerAMGSetMaxLevels(solver, 20, ierr) 250! conv. tolerance 251 call HYPRE_BoomerAMGSetTol(solver, 1.0d-7, ierr) 252 253! Now setup and solve! 254 call HYPRE_BoomerAMGSetup( 255 1 solver, parcsr_A, par_b, par_x, ierr) 256 call HYPRE_BoomerAMGSolve( 257 1 solver, parcsr_A, par_b, par_x, ierr) 258 259 260! Run info - needed logging turned on 261 call HYPRE_BoomerAMGGetNumIterations(solver, num_iterations, 262 1 ierr) 263 call HYPRE_BoomerAMGGetFinalReltvRes(solver, final_res_norm, 264 1 ierr) 265 266 267 if ( myid .eq. 0 ) then 268 print * 269 print '(A,I2)', " Iterations = ", num_iterations 270 print '(A,ES16.8)', 271 1 " Final Relative Residual Norm = ", final_res_norm 272 print * 273 endif 274 275! Destroy solver 276 call HYPRE_BoomerAMGDestroy(solver, ierr) 277 278! PCG (with DS) 279 elseif ( solver_id .eq. 50 ) then 280 281 282! Create solver 283 call HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver, ierr) 284 285! Set some parameters (See Reference Manual for more parameters) 286 call HYPRE_ParCSRPCGSetMaxIter(solver, 1000, ierr) 287 call HYPRE_ParCSRPCGSetTol(solver, 1.0d-7, ierr) 288 call HYPRE_ParCSRPCGSetTwoNorm(solver, 1, ierr) 289 call HYPRE_ParCSRPCGSetPrintLevel(solver, 2, ierr) 290 call HYPRE_ParCSRPCGSetLogging(solver, 1, ierr) 291 292! set ds (diagonal scaling) as the pcg preconditioner 293 precond_id = 1 294 call HYPRE_ParCSRPCGSetPrecond(solver, precond_id, 295 1 precond, ierr) 296 297 298 299! Now setup and solve! 300 call HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, 301 & par_x, ierr) 302 call HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, 303 & par_x, ierr) 304 305 306! Run info - needed logging turned on 307 308 call HYPRE_ParCSRPCGGetNumIterations(solver, num_iterations, 309 & ierr) 310 call HYPRE_ParCSRPCGGetFinalRelative(solver, final_res_norm, 311 & ierr) 312 313 if ( myid .eq. 0 ) then 314 print * 315 print *, "Iterations = ", num_iterations 316 print *, "Final Relative Residual Norm = ", final_res_norm 317 print * 318 endif 319 320! Destroy solver 321 call HYPRE_ParCSRPCGDestroy(solver, ierr) 322 323 324! PCG with AMG preconditioner 325 elseif ( solver_id == 1 ) then 326 327! Create solver 328 call HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver, ierr) 329 330! Set some parameters (See Reference Manual for more parameters) 331 call HYPRE_ParCSRPCGSetMaxIter(solver, 1000, ierr) 332 call HYPRE_ParCSRPCGSetTol(solver, 1.0d-7, ierr) 333 call HYPRE_ParCSRPCGSetTwoNorm(solver, 1, ierr) 334 call HYPRE_ParCSRPCGSetPrintLevel(solver, 2, ierr) 335 call HYPRE_ParCSRPCGSetLogging(solver, 1, ierr) 336 337! Now set up the AMG preconditioner and specify any parameters 338 339 call HYPRE_BoomerAMGCreate(precond, ierr) 340 341 342! Set some parameters (See Reference Manual for more parameters) 343 344! print less solver info since a preconditioner 345 call HYPRE_BoomerAMGSetPrintLevel(precond, 1, ierr); 346! Falgout coarsening 347 call HYPRE_BoomerAMGSetCoarsenType(precond, 6, ierr) 348! old defaults 349 call HYPRE_BoomerAMGSetOldDefault(precond, ierr) 350! SYMMETRIC G-S/Jacobi hybrid relaxation 351 call HYPRE_BoomerAMGSetRelaxType(precond, 6, ierr) 352! Sweeeps on each level 353 call HYPRE_BoomerAMGSetNumSweeps(precond, 1, ierr) 354! conv. tolerance 355 call HYPRE_BoomerAMGSetTol(precond, 0.0d0, ierr) 356! do only one iteration! 357 call HYPRE_BoomerAMGSetMaxIter(precond, 1, ierr) 358 359! set amg as the pcg preconditioner 360 precond_id = 2 361 call HYPRE_ParCSRPCGSetPrecond(solver, precond_id, 362 1 precond, ierr) 363 364 365! Now setup and solve! 366 call HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, 367 1 par_x, ierr) 368 call HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, 369 1 par_x, ierr) 370 371 372! Run info - needed logging turned on 373 374 call HYPRE_ParCSRPCGGetNumIterations(solver, num_iterations, 375 1 ierr) 376 call HYPRE_ParCSRPCGGetFinalRelative(solver, final_res_norm, 377 1 ierr) 378 379 if ( myid .eq. 0 ) then 380 print * 381 print *, "Iterations = ", num_iterations 382 print *, "Final Relative Residual Norm = ", final_res_norm 383 print * 384 endif 385 386! Destroy precond and solver 387 388 call HYPRE_BoomerAMGDestroy(precond, ierr) 389 call HYPRE_ParCSRPCGDestroy(solver, ierr) 390 391! PCG with ParaSails 392 elseif ( solver_id .eq. 8 ) then 393 394! Create solver 395 call HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, solver, ierr) 396 397! Set some parameters (See Reference Manual for more parameters) 398 call HYPRE_ParCSRPCGSetMaxIter(solver, 1000, ierr) 399 call HYPRE_ParCSRPCGSetTol(solver, 1.0d-7, ierr) 400 call HYPRE_ParCSRPCGSetTwoNorm(solver, 1, ierr) 401 call HYPRE_ParCSRPCGSetPrintLevel(solver, 2, ierr) 402 call HYPRE_ParCSRPCGSetLogging(solver, 1, ierr) 403 404! Now set up the Parasails preconditioner and specify any parameters 405 call HYPRE_ParaSailsCreate(MPI_COMM_WORLD, precond,ierr) 406 call HYPRE_ParaSailsSetParams(precond, 0.1d0, 1, ierr) 407 call HYPRE_ParaSailsSetFilter(precond, 0.05d0, ierr) 408 call HYPRE_ParaSailsSetSym(precond, 1, ierr) 409 call HYPRE_ParaSailsSetLogging(precond, 3, ierr) 410 411! set parsails as the pcg preconditioner 412 precond_id = 4 413 call HYPRE_ParCSRPCGSetPrecond(solver, precond_id, 414 1 precond, ierr) 415 416 417! Now setup and solve! 418 call HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, 419 1 par_x, ierr) 420 call HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, 421 1 par_x, ierr) 422 423 424! Run info - needed logging turned on 425 426 call HYPRE_ParCSRPCGGetNumIterations(solver, num_iterations, 427 1 ierr) 428 call HYPRE_ParCSRPCGGetFinalRelative(solver, final_res_norm, 429 1 ierr) 430 431 if ( myid .eq. 0 ) then 432 print * 433 print *, "Iterations = ", num_iterations 434 print *, "Final Relative Residual Norm = ", final_res_norm 435 print * 436 endif 437 438! Destroy precond and solver 439 440 call HYPRE_ParaSailsDestroy(precond, ierr) 441 call HYPRE_ParCSRPCGDestroy(solver, ierr) 442 443 else 444 if ( myid .eq. 0 ) then 445 print *,'Invalid solver id specified' 446 stop 447 endif 448 endif 449 450 451 452! Print the solution 453 if ( print_solution .ne. 0 ) then 454 call HYPRE_IJVectorPrint(x, "ij.out.x", ierr) 455 endif 456 457! Clean up 458 459 call HYPRE_IJMatrixDestroy(A, ierr) 460 call HYPRE_IJVectorDestroy(b, ierr) 461 call HYPRE_IJVectorDestroy(x, ierr) 462 463 call HYPRE_Finalize(ierr) 464 465! Finalize MPI 466 call MPI_Finalize(ierr) 467 468 stop 469 end 470