1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /* driver_commpkg.c*/
9 /* AHB 06/04 */
10 /* purpose: to test a new communication package for the ij interface */
11
12 /* 11/06 - if you want to use this, the the hypre_NewCommPkgCreate has to be
13 reinstated in parcsr_mv/new_commpkg.c - currently it won't compile*/
14
15
16 #include <stdlib.h>
17 #include <stdio.h>
18 #include <math.h>
19 /* #include <mpi.h> */
20
21 #include "_hypre_utilities.h"
22 #include "_hypre_parcsr_mv.h"
23 #include "HYPRE_parcsr_ls.h"
24
25 /* #include "_hypre_parcsr_ls.h"
26 #include "HYPRE.h"
27 #include "HYPRE_parcsr_mv.h"
28 #include "HYPRE_krylov.h" */
29
30
31
32 /*some debugging tools*/
33 #define mydebug 0
34 #define mpip_on 0
35
36 /*time an allgather in addition to the current commpkg -
37 since the allgather happens outside of the communication package.*/
38 #define time_gather 1
39
40 /* for timing multiple commpkg setup (if you want the time to be larger in the
41 hopes of getting smaller stds - often not effective) */
42 #define LOOP2 1
43
44
45 HYPRE_Int myBuildParLaplacian (HYPRE_Int argc , char *argv [], HYPRE_Int arg_index , HYPRE_ParCSRMatrix *A_ptr, HYPRE_Int parmprint );
46 HYPRE_Int myBuildParLaplacian27pt (HYPRE_Int argc , char *argv [], HYPRE_Int arg_index , HYPRE_ParCSRMatrix *A_ptr, HYPRE_Int parmprint );
47
48
49 void stats_mo(HYPRE_Real*, HYPRE_Int, HYPRE_Real *,HYPRE_Real *);
50
51 /*==========================================================================*/
52
53
54 /*------------------------------------------------------------------
55 *
56 * This tests an alternate comm package for ij
57 *
58 * options:
59 * -laplacian 3D 7pt stencil
60 * -27pt 3D 27pt laplacian
61 * -fromonecsrfile read matrix from a csr file
62 * -commpkg <HYPRE_Int> 1 = new comm. package
63 * 2 =old
64 * 3 = both (default)
65 * -loop <HYPRE_Int> number of times to loop (default is 0)
66 * -verbose print more error checking
67 * -noparmprint don't print the parameters
68 *-------------------------------------------------------------------*/
69
70
71 HYPRE_Int
main(HYPRE_Int argc,char * argv[])72 main( HYPRE_Int argc,
73 char *argv[] )
74 {
75
76
77 HYPRE_Int num_procs, myid;
78 HYPRE_Int verbose = 0, build_matrix_type = 1;
79 HYPRE_Int index, matrix_arg_index, commpkg_flag=3;
80 HYPRE_Int i, k, ierr=0;
81 HYPRE_Int row_start, row_end;
82 HYPRE_Int col_start, col_end, global_num_rows;
83 HYPRE_Int *row_part, *col_part;
84 char *csrfilename;
85 HYPRE_Int preload = 0, loop = 0, loop2 = LOOP2;
86 HYPRE_Int bcast_rows[2], *info;
87
88
89
90 hypre_ParCSRMatrix *parcsr_A, *small_A;
91 HYPRE_ParCSRMatrix A_temp, A_temp_small;
92 hypre_CSRMatrix *A_CSR;
93 hypre_ParCSRCommPkg *comm_pkg;
94
95
96 HYPRE_Int nx, ny, nz;
97 HYPRE_Int P, Q, R;
98 HYPRE_Int p, q, r;
99 HYPRE_Real values[4];
100
101 hypre_ParVector *x_new;
102 hypre_ParVector *y_new, *y;
103 HYPRE_Int *row_starts;
104 HYPRE_Real ans;
105 HYPRE_Real start_time, end_time, total_time, *loop_times;
106 HYPRE_Real T_avg, T_std;
107
108 HYPRE_Int noparmprint = 0;
109
110 #if mydebug
111 HYPRE_Int j, tmp_int;
112 #endif
113
114 /*-----------------------------------------------------------
115 * Initialize MPI
116 *-----------------------------------------------------------*/
117
118
119 hypre_MPI_Init(&argc, &argv);
120
121 hypre_MPI_Comm_size(hypre_MPI_COMM_WORLD, &num_procs );
122 hypre_MPI_Comm_rank(hypre_MPI_COMM_WORLD, &myid );
123
124
125
126 /*-----------------------------------------------------------
127 * default - is 27pt laplace
128 *-----------------------------------------------------------*/
129
130
131 build_matrix_type = 2;
132 matrix_arg_index = argc;
133
134 /*-----------------------------------------------------------
135 * Parse command line
136 *-----------------------------------------------------------*/
137
138 index = 1;
139 while ( index < argc)
140 {
141 if ( strcmp(argv[index], "-verbose") == 0 )
142 {
143 index++;
144 verbose = 1;
145 }
146 else if ( strcmp(argv[index], "-fromonecsrfile") == 0 )
147 {
148 index++;
149 build_matrix_type = 1;
150 matrix_arg_index = index; /*this tells where the name is*/
151 }
152 else if ( strcmp(argv[index], "-commpkg") == 0 )
153 {
154 index++;
155 commpkg_flag = atoi(argv[index++]);
156 }
157 else if ( strcmp(argv[index], "-laplacian") == 0 )
158 {
159 index++;
160 build_matrix_type = 2;
161 matrix_arg_index = index;
162 }
163 else if ( strcmp(argv[index], "-27pt") == 0 )
164 {
165 index++;
166 build_matrix_type = 4;
167 matrix_arg_index = index;
168 }
169 /*
170 else if ( strcmp(argv[index], "-nopreload") == 0 )
171 {
172 index++;
173 preload = 0;
174 }
175 */
176 else if ( strcmp(argv[index], "-loop") == 0 )
177 {
178 index++;
179 loop = atoi(argv[index++]);
180 }
181 else if ( strcmp(argv[index], "-noparmprint") == 0 )
182 {
183 index++;
184 noparmprint = 1;
185
186 }
187 else
188 {
189 index++;
190 /*hypre_printf("Warning: Unrecogized option '%s'\n",argv[index++] );*/
191 }
192 }
193
194
195
196 /*-----------------------------------------------------------
197 * Setup the Matrix problem
198 *-----------------------------------------------------------*/
199
200 /*-----------------------------------------------------------
201 * Get actual partitioning-
202 * read in an actual csr matrix.
203 *-----------------------------------------------------------*/
204
205
206 if (build_matrix_type ==1) /*read in a csr matrix from one file */
207 {
208 if (matrix_arg_index < argc)
209 {
210 csrfilename = argv[matrix_arg_index];
211 }
212 else
213 {
214 hypre_printf("Error: No filename specified \n");
215 exit(1);
216 }
217 if (myid == 0)
218 {
219 /*hypre_printf(" FromFile: %s\n", csrfilename);*/
220 A_CSR = hypre_CSRMatrixRead(csrfilename);
221 }
222 row_part = NULL;
223 col_part = NULL;
224
225 parcsr_A = hypre_CSRMatrixToParCSRMatrix(hypre_MPI_COMM_WORLD, A_CSR,
226 row_part, col_part);
227
228 if (myid == 0) hypre_CSRMatrixDestroy(A_CSR);
229 }
230 else if (build_matrix_type ==2)
231 {
232
233 myBuildParLaplacian(argc, argv, matrix_arg_index, &A_temp, !noparmprint);
234 parcsr_A = (hypre_ParCSRMatrix *) A_temp;
235
236 }
237 else if (build_matrix_type ==4)
238 {
239 myBuildParLaplacian27pt(argc, argv, matrix_arg_index, &A_temp, !noparmprint);
240 parcsr_A = (hypre_ParCSRMatrix *) A_temp;
241 }
242
243
244 /*-----------------------------------------------------------
245 * create a small problem so that timings are more accurate -
246 * code gets run twice (small laplace)
247 *-----------------------------------------------------------*/
248
249 /*this is no longer being used - preload = 0 is set at the beginning */
250
251 if (preload == 1)
252 {
253
254 /*hypre_printf("preload!\n");*/
255
256
257 values[1] = -1;
258 values[2] = -1;
259 values[3] = -1;
260 values[0] = - 6.0 ;
261
262 nx = 2;
263 ny = num_procs;
264 nz = 2;
265
266 P = 1;
267 Q = num_procs;
268 R = 1;
269
270 p = myid % P;
271 q = (( myid - p)/P) % Q;
272 r = ( myid - p - P*q)/( P*Q );
273
274 A_temp_small = (HYPRE_ParCSRMatrix) GenerateLaplacian(hypre_MPI_COMM_WORLD, nx, ny, nz,
275 P, Q, R, p, q, r, values);
276 small_A = (hypre_ParCSRMatrix *) A_temp_small;
277
278 /*do comm packages*/
279 hypre_NewCommPkgCreate(small_A);
280 hypre_NewCommPkgDestroy(small_A);
281
282 hypre_MatvecCommPkgCreate(small_A);
283 hypre_ParCSRMatrixDestroy(small_A);
284
285 }
286
287
288
289
290
291 /*-----------------------------------------------------------
292 * Prepare for timing
293 *-----------------------------------------------------------*/
294
295 /* instead of preloading, let's not time the first one if more than one*/
296
297
298 if (!loop)
299 {
300 loop = 1;
301 /* and don't do any timings */
302
303 }
304 else
305 {
306
307 loop +=1;
308 if (loop < 2) loop = 2;
309 }
310
311
312 loop_times = hypre_CTAlloc(HYPRE_Real, loop, HYPRE_MEMORY_HOST);
313
314
315
316 /******************************************************************************************/
317
318 hypre_MPI_Barrier(hypre_MPI_COMM_WORLD);
319
320 if (commpkg_flag == 1 || commpkg_flag ==3 )
321 {
322
323 /*-----------------------------------------------------------
324 * Create new comm package
325 *-----------------------------------------------------------*/
326
327
328
329 if (!myid) hypre_printf("********************************************************\n" );
330
331 /*do loop times*/
332 for (i=0; i< loop; i++)
333 {
334 loop_times[i] = 0.0;
335 for (k=0; k< loop2; k++)
336 {
337
338 hypre_MPI_Barrier(hypre_MPI_COMM_WORLD);
339
340 start_time = hypre_MPI_Wtime();
341
342 #if mpip_on
343 if (i==(loop-1)) hypre_MPI_Pcontrol(1);
344 #endif
345
346 hypre_NewCommPkgCreate(parcsr_A);
347
348 #if mpip_on
349 if (i==(loop-1)) hypre_MPI_Pcontrol(0);
350 #endif
351
352 end_time = hypre_MPI_Wtime();
353
354 end_time = end_time - start_time;
355
356 hypre_MPI_Allreduce(&end_time, &total_time, 1,
357 HYPRE_MPI_REAL, hypre_MPI_MAX, hypre_MPI_COMM_WORLD);
358
359 loop_times[i] += total_time;
360
361 if ( !((i+1)== loop && (k+1) == loop2)) hypre_NewCommPkgDestroy(parcsr_A);
362
363 }/*end of loop2 */
364
365
366 } /*end of loop*/
367
368
369
370 /* calculate the avg and std. */
371 if (loop > 1)
372 {
373
374 /* calculate the avg and std. */
375 stats_mo(loop_times, loop, &T_avg, &T_std);
376
377 if (!myid) hypre_printf(" NewCommPkgCreate: AVG. wall clock time = %f seconds\n", T_avg);
378 if (!myid) hypre_printf(" STD. for %d runs = %f\n", loop-1, T_std);
379 if (!myid) hypre_printf(" (Note: avg./std. timings exclude run 0.)\n");
380 if (!myid) hypre_printf("********************************************************\n" );
381 for (i=0; i< loop; i++)
382 {
383 if (!myid) hypre_printf(" run %d = %f sec.\n", i, loop_times[i]);
384 }
385 if (!myid) hypre_printf("********************************************************\n" );
386
387 }
388 else
389 {
390 if (!myid) hypre_printf("********************************************************\n" );
391 if (!myid) hypre_printf(" NewCommPkgCreate:\n");
392 if (!myid) hypre_printf(" run time = %f sec.\n", loop_times[0]);
393 if (!myid) hypre_printf("********************************************************\n" );
394 }
395
396
397 /*-----------------------------------------------------------
398 * Verbose printing
399 *-----------------------------------------------------------*/
400
401 /*some verification*/
402
403 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(parcsr_A);
404
405 if (verbose)
406 {
407
408 ierr = hypre_ParCSRMatrixGetLocalRange( parcsr_A,
409 &row_start, &row_end ,
410 &col_start, &col_end );
411
412
413 comm_pkg = hypre_ParCSRMatrixCommPkg(parcsr_A);
414
415 hypre_printf("myid = %i, my ACTUAL local range: [%i, %i]\n", myid,
416 row_start, row_end);
417
418
419 ierr = hypre_GetAssumedPartitionRowRange( myid, global_num_rows, &row_start,
420 &row_end);
421
422
423 hypre_printf("myid = %i, my assumed local range: [%i, %i]\n", myid,
424 row_start, row_end);
425
426 hypre_printf("myid = %d, num_recvs = %d\n", myid,
427 hypre_ParCSRCommPkgNumRecvs(comm_pkg) );
428
429 #if mydebug
430 for (i=0; i < hypre_ParCSRCommPkgNumRecvs(comm_pkg); i++)
431 {
432 hypre_printf("myid = %d, recv proc = %d, vec_starts = [%d : %d]\n",
433 myid, hypre_ParCSRCommPkgRecvProcs(comm_pkg)[i],
434 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg)[i],
435 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg)[i+1]-1);
436 }
437 #endif
438 hypre_printf("myid = %d, num_sends = %d\n", myid,
439 hypre_ParCSRCommPkgNumSends(comm_pkg) );
440
441 #if mydebug
442 for (i=0; i <hypre_ParCSRCommPkgNumSends(comm_pkg) ; i++)
443 {
444 tmp_int = hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i+1] -
445 hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i];
446 index = hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i];
447 for (j=0; j< tmp_int; j++)
448 {
449 hypre_printf("myid = %d, send proc = %d, send element = %d\n",myid,
450 hypre_ParCSRCommPkgSendProcs(comm_pkg)[i],
451 hypre_ParCSRCommPkgSendMapElmts(comm_pkg)[index+j]);
452 }
453 }
454 #endif
455 }
456 /*-----------------------------------------------------------
457 * To verify correctness (if commpkg_flag = 3)
458 *-----------------------------------------------------------*/
459
460 if (commpkg_flag == 3 )
461 {
462 /*do a matvec - we are assuming a square matrix */
463 row_starts = hypre_ParCSRMatrixRowStarts(parcsr_A);
464
465 x_new = hypre_ParVectorCreate(hypre_MPI_COMM_WORLD, global_num_rows, row_starts);
466 hypre_ParVectorInitialize(x_new);
467 hypre_ParVectorSetRandomValues(x_new, 1);
468
469 y_new = hypre_ParVectorCreate(hypre_MPI_COMM_WORLD, global_num_rows, row_starts);
470 hypre_ParVectorInitialize(y_new);
471 hypre_ParVectorSetConstantValues(y_new, 0.0);
472
473 /*y = 1.0*A*x+1.0*y */
474 hypre_ParCSRMatrixMatvec (1.0, parcsr_A, x_new, 1.0, y_new);
475 }
476
477 /*-----------------------------------------------------------
478 * Clean up after MyComm
479 *-----------------------------------------------------------*/
480
481
482 hypre_NewCommPkgDestroy(parcsr_A);
483
484 }
485
486
487
488
489
490
491 /******************************************************************************************/
492 /******************************************************************************************/
493
494 hypre_MPI_Barrier(hypre_MPI_COMM_WORLD);
495
496
497 if (commpkg_flag > 1 )
498 {
499
500 /*-----------------------------------------------------------
501 * Set up standard comm package
502 *-----------------------------------------------------------*/
503
504 bcast_rows[0] = 23;
505 bcast_rows[1] = 1789;
506
507 if (!myid) hypre_printf("********************************************************\n" );
508 /*do loop times*/
509 for (i=0; i< loop; i++)
510 {
511
512 loop_times[i] = 0.0;
513 for (k=0; k< loop2; k++)
514 {
515
516
517 hypre_MPI_Barrier(hypre_MPI_COMM_WORLD);
518
519
520 start_time = hypre_MPI_Wtime();
521
522 #if time_gather
523
524 info = hypre_CTAlloc(HYPRE_Int, num_procs, HYPRE_MEMORY_HOST);
525
526 hypre_MPI_Allgather(bcast_rows, 1, HYPRE_MPI_INT, info, 1, HYPRE_MPI_INT, hypre_MPI_COMM_WORLD);
527
528 #endif
529
530 hypre_MatvecCommPkgCreate(parcsr_A);
531
532 end_time = hypre_MPI_Wtime();
533
534
535 end_time = end_time - start_time;
536
537 hypre_MPI_Allreduce(&end_time, &total_time, 1,
538 HYPRE_MPI_REAL, hypre_MPI_MAX, hypre_MPI_COMM_WORLD);
539
540 loop_times[i] += total_time;
541
542
543 if ( !((i+1)== loop && (k+1) == loop2)) hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(parcsr_A));
544
545 }/* end of loop 2*/
546
547
548 } /*end of loop*/
549
550 /* calculate the avg and std. */
551 if (loop > 1)
552 {
553
554 stats_mo(loop_times, loop, &T_avg, &T_std);
555 if (!myid) hypre_printf("Current CommPkgCreate: AVG. wall clock time = %f seconds\n", T_avg);
556 if (!myid) hypre_printf(" STD. for %d runs = %f\n", loop-1, T_std);
557 if (!myid) hypre_printf(" (Note: avg./std. timings exclude run 0.)\n");
558 if (!myid) hypre_printf("********************************************************\n" );
559 for (i=0; i< loop; i++)
560 {
561 if (!myid) hypre_printf(" run %d = %f sec.\n", i, loop_times[i]);
562 }
563 if (!myid) hypre_printf("********************************************************\n" );
564
565 }
566 else
567 {
568 if (!myid) hypre_printf("********************************************************\n" );
569 if (!myid) hypre_printf(" Current CommPkgCreate:\n");
570 if (!myid) hypre_printf(" run time = %f sec.\n", loop_times[0]);
571 if (!myid) hypre_printf("********************************************************\n" );
572 }
573
574
575
576
577
578 /*-----------------------------------------------------------
579 * Verbose printing
580 *-----------------------------------------------------------*/
581
582 /*some verification*/
583
584
585 if (verbose)
586 {
587
588 ierr = hypre_ParCSRMatrixGetLocalRange( parcsr_A,
589 &row_start, &row_end ,
590 &col_start, &col_end );
591
592
593 comm_pkg = hypre_ParCSRMatrixCommPkg(parcsr_A);
594
595 hypre_printf("myid = %i, std - my local range: [%i, %i]\n", myid,
596 row_start, row_end);
597
598 ierr = hypre_ParCSRMatrixGetLocalRange( parcsr_A,
599 &row_start, &row_end ,
600 &col_start, &col_end );
601
602 hypre_printf("myid = %d, std - num_recvs = %d\n", myid,
603 hypre_ParCSRCommPkgNumRecvs(comm_pkg) );
604
605 #if mydebug
606 for (i=0; i < hypre_ParCSRCommPkgNumRecvs(comm_pkg); i++)
607 {
608 hypre_printf("myid = %d, std - recv proc = %d, vec_starts = [%d : %d]\n",
609 myid, hypre_ParCSRCommPkgRecvProcs(comm_pkg)[i],
610 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg)[i],
611 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg)[i+1]-1);
612 }
613 #endif
614 hypre_printf("myid = %d, std - num_sends = %d\n", myid,
615 hypre_ParCSRCommPkgNumSends(comm_pkg));
616
617
618 #if mydebug
619 for (i=0; i <hypre_ParCSRCommPkgNumSends(comm_pkg) ; i++)
620 {
621 tmp_int = hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i+1] -
622 hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i];
623 index = hypre_ParCSRCommPkgSendMapStarts(comm_pkg)[i];
624 for (j=0; j< tmp_int; j++)
625 {
626 hypre_printf("myid = %d, std - send proc = %d, send element = %d\n",myid,
627 hypre_ParCSRCommPkgSendProcs(comm_pkg)[i],
628 hypre_ParCSRCommPkgSendMapElmts(comm_pkg)[index+j]);
629 }
630 }
631 #endif
632 }
633
634 /*-----------------------------------------------------------
635 * Verify correctness
636 *-----------------------------------------------------------*/
637
638
639
640 if (commpkg_flag == 3 )
641 {
642 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(parcsr_A);
643 row_starts = hypre_ParCSRMatrixRowStarts(parcsr_A);
644
645
646 y = hypre_ParVectorCreate(hypre_MPI_COMM_WORLD, global_num_rows,row_starts);
647 hypre_ParVectorInitialize(y);
648 hypre_ParVectorSetConstantValues(y, 0.0);
649
650 hypre_ParCSRMatrixMatvec (1.0, parcsr_A, x_new, 1.0, y);
651 }
652 }
653
654
655
656
657
658
659 /*-----------------------------------------------------------
660 * Compare matvecs for both comm packages (3)
661 *-----------------------------------------------------------*/
662
663 if (commpkg_flag == 3 )
664 {
665 /*make sure that y and y_new are the same - now y_new should=0*/
666 hypre_ParVectorAxpy( -1.0, y, y_new );
667
668
669 hypre_ParVectorSetRandomValues(y, 1);
670
671 ans = hypre_ParVectorInnerProd( y, y_new );
672 if (!myid)
673 {
674
675 if ( fabs(ans) > 1e-8 )
676 {
677 hypre_printf("!!!!! WARNING !!!!! should be zero if correct = %6.10f\n",
678 ans);
679 }
680 else
681 {
682 hypre_printf("Matvecs match ( should be zero = %6.10f )\n",
683 ans);
684 }
685 }
686
687
688 }
689
690
691 /*-----------------------------------------------------------
692 * Clean up
693 *-----------------------------------------------------------*/
694
695
696 hypre_ParCSRMatrixDestroy(parcsr_A); /*this calls the standard comm
697 package destroy - but we'll destroy
698 ours separately until it is
699 incorporated */
700
701 if (commpkg_flag == 3 )
702 {
703
704 hypre_ParVectorDestroy(x_new);
705 hypre_ParVectorDestroy(y);
706 hypre_ParVectorDestroy(y_new);
707 }
708
709
710
711
712 hypre_MPI_Finalize();
713
714 return(ierr);
715
716
717 }
718
719
720
721
722
723 /*------------------------------------
724 * Calculate the average and STD
725 * throw away 1st timing
726 *------------------------------------*/
727
stats_mo(HYPRE_Real array[],HYPRE_Int n,HYPRE_Real * Tavg,HYPRE_Real * Tstd)728 void stats_mo(HYPRE_Real array[], HYPRE_Int n, HYPRE_Real *Tavg,HYPRE_Real *Tstd)
729 {
730
731 HYPRE_Int i;
732 HYPRE_Real atmp, tmp=0.0;
733 HYPRE_Real avg = 0.0, std;
734
735
736 for(i=1; i<n; i++) {
737 atmp = array[i];
738 avg += atmp;
739 tmp += atmp*atmp;
740 }
741
742 n = n-1;
743 avg = avg/(HYPRE_Real) n;
744 tmp = tmp/(HYPRE_Real) n;
745
746 tmp = fabs(tmp - avg*avg);
747 std = sqrt(tmp);
748
749 *Tavg = avg;
750 *Tstd = std;
751 }
752
753
754
755 /*These next two functions are from ij.c in linear_solvers/tests */
756
757
758 /*----------------------------------------------------------------------
759 * Build 27-point laplacian in 3D,
760 * Parameters given in command line.
761 *----------------------------------------------------------------------*/
762
763 HYPRE_Int
myBuildParLaplacian27pt(HYPRE_Int argc,char * argv[],HYPRE_Int arg_index,HYPRE_ParCSRMatrix * A_ptr,HYPRE_Int parmprint)764 myBuildParLaplacian27pt( HYPRE_Int argc,
765 char *argv[],
766 HYPRE_Int arg_index,
767 HYPRE_ParCSRMatrix *A_ptr , HYPRE_Int parmprint )
768 {
769 HYPRE_Int nx, ny, nz;
770 HYPRE_Int P, Q, R;
771
772 HYPRE_ParCSRMatrix A;
773
774 HYPRE_Int num_procs, myid;
775 HYPRE_Int p, q, r;
776 HYPRE_Real *values;
777
778 /*-----------------------------------------------------------
779 * Initialize some stuff
780 *-----------------------------------------------------------*/
781
782 hypre_MPI_Comm_size(hypre_MPI_COMM_WORLD, &num_procs );
783 hypre_MPI_Comm_rank(hypre_MPI_COMM_WORLD, &myid );
784
785 /*-----------------------------------------------------------
786 * Set defaults
787 *-----------------------------------------------------------*/
788
789 nx = 10;
790 ny = 10;
791 nz = 10;
792
793 P = 1;
794 Q = num_procs;
795 R = 1;
796
797 /*-----------------------------------------------------------
798 * Parse command line
799 *-----------------------------------------------------------*/
800 arg_index = 0;
801 while (arg_index < argc)
802 {
803 if ( strcmp(argv[arg_index], "-n") == 0 )
804 {
805 arg_index++;
806 nx = atoi(argv[arg_index++]);
807 ny = atoi(argv[arg_index++]);
808 nz = atoi(argv[arg_index++]);
809 }
810 else if ( strcmp(argv[arg_index], "-P") == 0 )
811 {
812 arg_index++;
813 P = atoi(argv[arg_index++]);
814 Q = atoi(argv[arg_index++]);
815 R = atoi(argv[arg_index++]);
816 }
817 else
818 {
819 arg_index++;
820 }
821 }
822
823 /*-----------------------------------------------------------
824 * Check a few things
825 *-----------------------------------------------------------*/
826
827 if ((P*Q*R) != num_procs)
828 {
829 hypre_printf("Error: Invalid number of processors or processor topology \n");
830 exit(1);
831 }
832
833 /*-----------------------------------------------------------
834 * Print driver parameters
835 *-----------------------------------------------------------*/
836
837 if (myid == 0 && parmprint)
838 {
839 hypre_printf(" Laplacian_27pt:\n");
840 hypre_printf(" (nx, ny, nz) = (%d, %d, %d)\n", nx, ny, nz);
841 hypre_printf(" (Px, Py, Pz) = (%d, %d, %d)\n\n", P, Q, R);
842 }
843
844 /*-----------------------------------------------------------
845 * Set up the grid structure
846 *-----------------------------------------------------------*/
847
848 /* compute p,q,r from P,Q,R and myid */
849 p = myid % P;
850 q = (( myid - p)/P) % Q;
851 r = ( myid - p - P*q)/( P*Q );
852
853 /*-----------------------------------------------------------
854 * Generate the matrix
855 *-----------------------------------------------------------*/
856
857 values = hypre_CTAlloc(HYPRE_Real, 2, HYPRE_MEMORY_HOST);
858
859 values[0] = 26.0;
860 if (nx == 1 || ny == 1 || nz == 1)
861 values[0] = 8.0;
862 if (nx*ny == 1 || nx*nz == 1 || ny*nz == 1)
863 values[0] = 2.0;
864 values[1] = -1.;
865
866 A = (HYPRE_ParCSRMatrix) GenerateLaplacian27pt(hypre_MPI_COMM_WORLD,
867 nx, ny, nz, P, Q, R, p, q, r, values);
868
869 hypre_TFree(values, HYPRE_MEMORY_HOST);
870
871 *A_ptr = A;
872
873 return (0);
874 }
875
876
877 /*----------------------------------------------------------------------
878 * Build standard 7-point laplacian in 3D with grid and anisotropy.
879 * Parameters given in command line.
880 *----------------------------------------------------------------------*/
881
882
883 HYPRE_Int
myBuildParLaplacian(HYPRE_Int argc,char * argv[],HYPRE_Int arg_index,HYPRE_ParCSRMatrix * A_ptr,HYPRE_Int parmprint)884 myBuildParLaplacian( HYPRE_Int argc,
885 char *argv[],
886 HYPRE_Int arg_index,
887 HYPRE_ParCSRMatrix *A_ptr , HYPRE_Int parmprint )
888 {
889 HYPRE_Int nx, ny, nz;
890 HYPRE_Int P, Q, R;
891 HYPRE_Real cx, cy, cz;
892
893 HYPRE_ParCSRMatrix A;
894
895 HYPRE_Int num_procs, myid;
896 HYPRE_Int p, q, r;
897 HYPRE_Real *values;
898
899 /*-----------------------------------------------------------
900 * Initialize some stuff
901 *-----------------------------------------------------------*/
902
903 hypre_MPI_Comm_size(hypre_MPI_COMM_WORLD, &num_procs );
904 hypre_MPI_Comm_rank(hypre_MPI_COMM_WORLD, &myid );
905
906 /*-----------------------------------------------------------
907 * Set defaults
908 *-----------------------------------------------------------*/
909
910 nx = 10;
911 ny = 10;
912 nz = 10;
913
914 P = 1;
915 Q = num_procs;
916 R = 1;
917
918 cx = 1.;
919 cy = 1.;
920 cz = 1.;
921
922 /*-----------------------------------------------------------
923 * Parse command line
924 *-----------------------------------------------------------*/
925 arg_index = 0;
926 while (arg_index < argc)
927 {
928 if ( strcmp(argv[arg_index], "-n") == 0 )
929 {
930 arg_index++;
931 nx = atoi(argv[arg_index++]);
932 ny = atoi(argv[arg_index++]);
933 nz = atoi(argv[arg_index++]);
934 }
935 else if ( strcmp(argv[arg_index], "-P") == 0 )
936 {
937 arg_index++;
938 P = atoi(argv[arg_index++]);
939 Q = atoi(argv[arg_index++]);
940 R = atoi(argv[arg_index++]);
941 }
942 else if ( strcmp(argv[arg_index], "-c") == 0 )
943 {
944 arg_index++;
945 cx = atof(argv[arg_index++]);
946 cy = atof(argv[arg_index++]);
947 cz = atof(argv[arg_index++]);
948 }
949 else
950 {
951 arg_index++;
952 }
953 }
954
955 /*-----------------------------------------------------------
956 * Check a few things
957 *-----------------------------------------------------------*/
958
959 if ((P*Q*R) != num_procs)
960 {
961 hypre_printf("Error: Invalid number of processors or processor topology \n");
962 exit(1);
963 }
964
965 /*-----------------------------------------------------------
966 * Print driver parameters
967 *-----------------------------------------------------------*/
968
969 if (myid == 0 && parmprint)
970 {
971 hypre_printf(" Laplacian:\n");
972 hypre_printf(" (nx, ny, nz) = (%d, %d, %d)\n", nx, ny, nz);
973 hypre_printf(" (Px, Py, Pz) = (%d, %d, %d)\n", P, Q, R);
974 hypre_printf(" (cx, cy, cz) = (%f, %f, %f)\n\n", cx, cy, cz);
975 }
976
977 /*-----------------------------------------------------------
978 * Set up the grid structure
979 *-----------------------------------------------------------*/
980
981 /* compute p,q,r from P,Q,R and myid */
982 p = myid % P;
983 q = (( myid - p)/P) % Q;
984 r = ( myid - p - P*q)/( P*Q );
985
986 /*-----------------------------------------------------------
987 * Generate the matrix
988 *-----------------------------------------------------------*/
989
990 values = hypre_CTAlloc(HYPRE_Real, 4, HYPRE_MEMORY_HOST);
991
992 values[1] = -cx;
993 values[2] = -cy;
994 values[3] = -cz;
995
996 values[0] = 0.;
997 if (nx > 1)
998 {
999 values[0] += 2.0*cx;
1000 }
1001 if (ny > 1)
1002 {
1003 values[0] += 2.0*cy;
1004 }
1005 if (nz > 1)
1006 {
1007 values[0] += 2.0*cz;
1008 }
1009
1010 A = (HYPRE_ParCSRMatrix) GenerateLaplacian(hypre_MPI_COMM_WORLD, nx, ny, nz,
1011 P, Q, R, p, q, r, values);
1012
1013 hypre_TFree(values, HYPRE_MEMORY_HOST);
1014
1015
1016 *A_ptr = A;
1017
1018 return (0);
1019 }
1020