12 #include "../config/config.hpp"
17 #include "../fem/fem.hpp"
29 HypreParVector::HypreParVector(MPI_Comm comm,
int glob_size,
32 x = hypre_ParVectorCreate(comm,glob_size,col);
33 hypre_ParVectorInitialize(x);
34 hypre_ParVectorSetPartitioningOwner(x,0);
36 hypre_ParVectorSetDataOwner(x,1);
37 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
39 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
44 double *_data,
int *col) :
Vector()
46 x = hypre_ParVectorCreate(comm,glob_size,col);
47 hypre_ParVectorSetDataOwner(x,1);
48 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),0);
49 hypre_ParVectorSetPartitioningOwner(x,0);
50 hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
53 hypre_ParVectorInitialize(x);
55 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
61 x = hypre_ParVectorCreate(y.x -> comm, y.x -> global_size,
63 hypre_ParVectorInitialize(x);
64 hypre_ParVectorSetPartitioningOwner(x,0);
65 hypre_ParVectorSetDataOwner(x,1);
66 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
68 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
75 x = hypre_ParVectorInDomainOf(A);
77 x = hypre_ParVectorInRangeOf(A);
79 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
85 x = (hypre_ParVector *) y;
87 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
95 hypre_ParVectorInitialize(x);
96 hypre_ParVectorSetPartitioningOwner(x,0);
98 hypre_ParVectorSetDataOwner(x,1);
99 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
101 hypre_VectorSize(hypre_ParVectorLocalVector(x)));
105 HypreParVector::operator hypre_ParVector*()
const
110 #ifndef HYPRE_PAR_VECTOR_STRUCT
111 HypreParVector::operator HYPRE_ParVector()
const
113 return (HYPRE_ParVector) x;
119 hypre_Vector *hv = hypre_ParVectorToVectorAll(*
this);
122 hypre_SeqVectorSetDataOwner(hv,0);
123 hypre_SeqVectorDestroy(hv);
129 hypre_ParVectorSetConstantValues(x,d);
140 for (
int i = 0; i <
size; i++)
147 Vector::data = hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
152 return hypre_ParVectorSetRandomValues(x,seed);
157 hypre_ParVectorPrint(x,fname);
163 hypre_ParVectorDestroy(x);
169 return hypre_ParVectorInnerProd(*x, *y);
174 return hypre_ParVectorInnerProd(x, y);
180 :
Operator(diag->Height(), diag->Width())
182 A = hypre_ParCSRMatrixCreate(comm, glob_size, glob_size, row_starts,
184 hypre_ParCSRMatrixSetDataOwner(A,0);
185 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
186 hypre_ParCSRMatrixSetColStartsOwner(A,0);
188 hypre_CSRMatrixSetDataOwner(A->diag,0);
189 hypre_CSRMatrixI(A->diag) = diag->
GetI();
190 hypre_CSRMatrixJ(A->diag) = diag->
GetJ();
191 hypre_CSRMatrixData(A->diag) = diag->
GetData();
192 hypre_CSRMatrixSetRownnz(A->diag);
194 hypre_CSRMatrixSetDataOwner(A->offd,1);
195 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(
int, diag->
Height()+1);
202 hypre_ParCSRMatrixSetNumNonzeros(A);
205 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
207 hypre_MatvecCommPkgCreate(A);
215 int global_num_rows,
int global_num_cols,
216 int *row_starts,
int *col_starts,
218 :
Operator(diag->Height(), diag->Width())
220 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
221 row_starts, col_starts,
223 hypre_ParCSRMatrixSetDataOwner(A,0);
224 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
225 hypre_ParCSRMatrixSetColStartsOwner(A,0);
227 hypre_CSRMatrixSetDataOwner(A->diag,0);
228 hypre_CSRMatrixI(A->diag) = diag->
GetI();
229 hypre_CSRMatrixJ(A->diag) = diag->
GetJ();
230 hypre_CSRMatrixData(A->diag) = diag->
GetData();
231 hypre_CSRMatrixSetRownnz(A->diag);
233 hypre_CSRMatrixSetDataOwner(A->offd,1);
234 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(
int, diag->
Height()+1);
236 hypre_ParCSRMatrixSetNumNonzeros(A);
239 if (row_starts == col_starts)
240 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
242 hypre_MatvecCommPkgCreate(A);
249 int global_num_rows,
int global_num_cols,
250 int *row_starts,
int *col_starts,
253 :
Operator(diag->Height(), diag->Width())
255 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
256 row_starts, col_starts,
259 hypre_ParCSRMatrixSetDataOwner(A,0);
260 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
261 hypre_ParCSRMatrixSetColStartsOwner(A,0);
263 hypre_CSRMatrixSetDataOwner(A->diag,0);
264 hypre_CSRMatrixI(A->diag) = diag->
GetI();
265 hypre_CSRMatrixJ(A->diag) = diag->
GetJ();
266 hypre_CSRMatrixData(A->diag) = diag->
GetData();
267 hypre_CSRMatrixSetRownnz(A->diag);
269 hypre_CSRMatrixSetDataOwner(A->offd,0);
270 hypre_CSRMatrixI(A->offd) = offd->
GetI();
271 hypre_CSRMatrixJ(A->offd) = offd->
GetJ();
272 hypre_CSRMatrixData(A->offd) = offd->
GetData();
273 hypre_CSRMatrixSetRownnz(A->offd);
275 hypre_ParCSRMatrixColMapOffd(A) = cmap;
277 hypre_ParCSRMatrixSetNumNonzeros(A);
280 if (row_starts == col_starts)
281 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
283 hypre_MatvecCommPkgCreate(A);
294 mfem_error(
"HypreParMatrix::HypreParMatrix: sm_a==NULL");
297 hypre_CSRMatrix *csr_a;
299 csr_a = hypre_CSRMatrixCreate(sm_a ->
Height(), sm_a ->
Width(),
300 sm_a -> NumNonZeroElems());
302 hypre_CSRMatrixSetDataOwner(csr_a,0);
303 hypre_CSRMatrixI(csr_a) = sm_a -> GetI();
304 hypre_CSRMatrixJ(csr_a) = sm_a -> GetJ();
305 hypre_CSRMatrixData(csr_a) = sm_a -> GetData();
306 hypre_CSRMatrixSetRownnz(csr_a);
308 A = hypre_CSRMatrixToParCSRMatrix(comm, csr_a, row_starts, col_starts);
317 if (row_starts == col_starts)
318 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
320 hypre_MatvecCommPkgCreate(A);
324 int global_num_rows,
int global_num_cols,
325 int *row_starts,
int *col_starts,
Table *diag)
328 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
329 row_starts, col_starts, 0, nnz, 0);
330 hypre_ParCSRMatrixSetDataOwner(A,1);
331 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
332 hypre_ParCSRMatrixSetColStartsOwner(A,0);
334 hypre_CSRMatrixSetDataOwner(A->diag,1);
335 hypre_CSRMatrixI(A->diag) = diag->
GetI();
336 hypre_CSRMatrixJ(A->diag) = diag->
GetJ();
338 hypre_CSRMatrixData(A->diag) = hypre_TAlloc(
double, nnz);
339 for (
int k = 0; k < nnz; k++)
340 (hypre_CSRMatrixData(A->diag))[k] = 1.0;
342 hypre_CSRMatrixSetRownnz(A->diag);
344 hypre_CSRMatrixSetDataOwner(A->offd,1);
345 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(
int, diag->
Size()+1);
347 hypre_ParCSRMatrixSetNumNonzeros(A);
350 if (row_starts == col_starts)
351 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
353 hypre_MatvecCommPkgCreate(A);
364 int *i_diag,
int *j_diag,
365 int *i_offd,
int *j_offd,
366 int *cmap,
int cmap_size)
368 int diag_col, offd_col;
370 if (HYPRE_AssumedPartitionCheck())
372 diag_col = i_diag[row[1]-row[0]];
373 offd_col = i_offd[row[1]-row[0]];
375 A = hypre_ParCSRMatrixCreate(comm, row[2], col[2], row, col,
376 cmap_size, diag_col, offd_col);
380 diag_col = i_diag[row[
id+1]-row[id]];
381 offd_col = i_offd[row[
id+1]-row[id]];
383 A = hypre_ParCSRMatrixCreate(comm, row[np], col[np], row, col,
384 cmap_size, diag_col, offd_col);
387 hypre_ParCSRMatrixSetDataOwner(A,1);
388 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
389 hypre_ParCSRMatrixSetColStartsOwner(A,0);
393 double *a_diag = hypre_TAlloc(
double, diag_col);
394 for (i = 0; i < diag_col; i++)
397 double *a_offd = hypre_TAlloc(
double, offd_col);
398 for (i = 0; i < offd_col; i++)
401 hypre_CSRMatrixSetDataOwner(A->diag,1);
402 hypre_CSRMatrixI(A->diag) = i_diag;
403 hypre_CSRMatrixJ(A->diag) = j_diag;
404 hypre_CSRMatrixData(A->diag) = a_diag;
405 hypre_CSRMatrixSetRownnz(A->diag);
407 hypre_CSRMatrixSetDataOwner(A->offd,1);
408 hypre_CSRMatrixI(A->offd) = i_offd;
409 hypre_CSRMatrixJ(A->offd) = j_offd;
410 hypre_CSRMatrixData(A->offd) = a_offd;
411 hypre_CSRMatrixSetRownnz(A->offd);
413 hypre_ParCSRMatrixColMapOffd(A) = cmap;
415 hypre_ParCSRMatrixSetNumNonzeros(A);
419 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
421 hypre_MatvecCommPkgCreate(A);
431 int glob_ncols,
int *I,
int *J,
double *data,
432 int *rows,
int *cols)
436 hypre_CSRMatrix *local = hypre_CSRMatrixCreate(nrows, glob_ncols, nnz);
437 hypre_CSRMatrixI(local) = I;
438 hypre_CSRMatrixJ(local) = J;
439 hypre_CSRMatrixData(local) = data;
440 hypre_CSRMatrixRownnz(local) = NULL;
441 hypre_CSRMatrixOwnsData(local) = 1;
442 hypre_CSRMatrixNumRownnz(local) = nrows;
445 if (HYPRE_AssumedPartitionCheck())
452 MPI_Comm_rank(comm, &myid);
453 MPI_Comm_size(comm, &part_size);
458 int *row_starts = hypre_TAlloc(
int, part_size);
459 int *col_starts = hypre_TAlloc(
int, part_size);
460 for (
int i = 0; i < part_size; i++)
462 row_starts[i] = rows[i];
463 col_starts[i] = cols[i];
467 A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
468 row_starts, col_starts, 0, 0, 0);
469 hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
470 hypre_ParCSRMatrixOwnsColStarts(A) = 1;
471 GenerateDiagAndOffd(local, A, col_starts[myid], col_starts[myid+1]-1);
472 hypre_ParCSRMatrixSetNumNonzeros(A);
475 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
476 hypre_MatvecCommPkgCreate(A);
479 hypre_CSRMatrixI(local) = NULL;
480 hypre_CSRMatrixJ(local) = NULL;
481 hypre_CSRMatrixData(local) = NULL;
482 hypre_CSRMatrixDestroy(local);
494 if (hypre_ParCSRMatrixCommPkg(A))
495 hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(A));
497 hypre_ParCSRMatrixCommPkg(A) = comm_pkg;
503 if (CommPkg == NULL || CommPkg != hypre_ParCSRMatrixCommPkg(A))
504 mfem_error(
"\nHypreParMatrix::CheckCommPkg()");
512 hypre_TFree(CommPkg->send_procs);
513 hypre_TFree(CommPkg->send_map_starts);
514 hypre_TFree(CommPkg->send_map_elmts);
515 hypre_TFree(CommPkg->recv_procs);
516 hypre_TFree(CommPkg->recv_vec_starts);
517 if (CommPkg->send_mpi_types)
518 hypre_TFree(CommPkg->send_mpi_types);
519 if (CommPkg->recv_mpi_types)
520 hypre_TFree(CommPkg->recv_mpi_types);
521 if (hypre_ParCSRMatrixCommPkg(A) == CommPkg)
522 hypre_ParCSRMatrixCommPkg(A) = NULL;
527 HypreParMatrix::operator hypre_ParCSRMatrix*()
529 return (
this) ? A : NULL;
532 #ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
533 HypreParMatrix::operator HYPRE_ParCSRMatrix()
535 return (
this) ? (HYPRE_ParCSRMatrix) A : (HYPRE_ParCSRMatrix) NULL;
541 hypre_ParCSRMatrix *R = A;
548 int size=hypre_CSRMatrixNumRows(A->diag);
550 for (
int j = 0; j < size; j++)
552 diag(j) = A->diag->data[A->diag->i[j]];
554 if (A->diag->j[A->diag->i[j]] != j)
563 hypre_ParCSRMatrix * At;
564 hypre_ParCSRMatrixTranspose(A, &At, 1);
565 hypre_ParCSRMatrixSetNumNonzeros(At);
567 hypre_MatvecCommPkgCreate(At);
575 return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
597 hypre_ParCSRMatrixMatvec(a, A, *X, b, *Y);
601 double b,
Vector &y)
const
622 hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
628 return hypre_ParCSRMatrixMatvec(a,A,(hypre_ParVector *)x,b,(hypre_ParVector *)y);
634 return hypre_ParCSRMatrixMatvecT(a,A,x,b,y);
640 if(hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
643 if(hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
644 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
646 int size=hypre_CSRMatrixNumRows(A->diag);
647 double *Adiag_data = hypre_CSRMatrixData(A->diag);
648 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
651 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
652 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
655 for (
int i(0); i < size; ++i)
658 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
659 Adiag_data[jj] *= val;
660 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
661 Aoffd_data[jj] *= val;
668 if(hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
671 if(hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
672 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
674 int size=hypre_CSRMatrixNumRows(A->diag);
675 double *Adiag_data = hypre_CSRMatrixData(A->diag);
676 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
679 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
680 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
683 for (
int i(0); i < size; ++i)
687 mfem_error(
"HypreParMatrix::InvDiagScale : Division by 0");
690 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
691 Adiag_data[jj] *= val;
692 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
693 Aoffd_data[jj] *= val;
699 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
702 int size=hypre_CSRMatrixNumRows(A->diag);
705 double *Adiag_data = hypre_CSRMatrixData(A->diag);
706 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
707 for (jj = 0; jj < Adiag_i[size]; ++jj)
710 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
711 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
712 for (jj = 0; jj < Aoffd_i[size]; ++jj)
718 hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
723 if (A) hypre_ParCSRMatrixDestroy(A);
725 hypre_ParCSRMatrixReadIJ(comm, fname, &io, &jo, &A);
726 hypre_ParCSRMatrixSetNumNonzeros(A);
728 hypre_MatvecCommPkgCreate(A);
737 if (hypre_ParCSRMatrixCommPkg(A))
738 hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(A));
739 hypre_ParCSRMatrixCommPkg(A) = NULL;
741 if (hypre_CSRMatrixOwnsData(A->diag))
743 hypre_CSRMatrixDestroy(A->diag);
748 if (hypre_CSRMatrixRownnz(A->diag))
749 hypre_TFree(hypre_CSRMatrixRownnz(A->diag));
750 hypre_TFree(A->diag);
754 if (hypre_CSRMatrixOwnsData(A->offd))
756 hypre_CSRMatrixDestroy(A->offd);
761 if (hypre_CSRMatrixRownnz(A->offd))
762 hypre_TFree(hypre_CSRMatrixRownnz(A->offd));
764 hypre_ParCSRMatrixDestroy(A);
773 hypre_ParCSRMatrix * ab;
774 ab = hypre_ParMatmul(*A,*B);
776 hypre_MatvecCommPkgCreate(ab);
783 int P_owns_its_col_starts =
784 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
785 hypre_ParCSRMatrix * rap;
786 hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
788 hypre_ParCSRMatrixSetNumNonzeros(rap);
790 if (!P_owns_its_col_starts)
794 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
795 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
802 int P_owns_its_col_starts =
803 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
804 int Rt_owns_its_col_starts =
805 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
807 hypre_ParCSRMatrix * rap;
808 hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
810 hypre_ParCSRMatrixSetNumNonzeros(rap);
812 if (!P_owns_its_col_starts)
816 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
818 if (!Rt_owns_its_col_starts)
822 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
832 Ae.
Mult(x, b, -1.0, 1.0);
834 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *)A);
835 double *data = hypre_CSRMatrixData(A_diag);
836 int *I = hypre_CSRMatrixI(A_diag);
838 int *J = hypre_CSRMatrixJ(A_diag);
840 hypre_CSRMatrixI(hypre_ParCSRMatrixOffd((hypre_ParCSRMatrix *)A));
843 for (
int i = 0; i < ess_dof_list.
Size(); i++)
845 int r = ess_dof_list[i];
846 b(r) = data[I[r]] * x(r);
850 if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
869 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
870 int num_rows = hypre_CSRMatrixNumRows(A_diag);
872 double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
873 double *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
875 for (
int i = 0; i < N; i++)
878 hypre_ParVectorCopy(f, r);
879 hypre_ParCSRMatrixMatvec(-1.0, A, u, 1.0, r);
882 (0 == (i % 2)) ? coef = lambda : coef = mu;
884 for (
int i = 0; i < num_rows; i++)
886 u_data[i] += coef*r_data[i] / max_eig;
908 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
909 int num_rows = hypre_CSRMatrixNumRows(A_diag);
911 double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
913 double *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
914 double *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
915 double *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
916 double *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
918 hypre_ParVectorCopy(u, x0);
921 hypre_ParVectorCopy(f, x1);
922 hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
924 for (
int i = 0; i < num_rows; i++)
926 x1_data[i] /= -max_eig;
930 for (
int i = 0; i < num_rows; i++)
932 x1_data[i] = x0_data[i] -x1_data[i];
936 for (
int i = 0; i < num_rows; i++)
938 x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
941 for (
int n = 2; n <= poly_order; n++)
944 hypre_ParVectorCopy(f, x2);
945 hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
947 for (
int i = 0; i < num_rows; i++)
949 x2_data[i] /= -max_eig;
957 for (
int i = 0; i < num_rows; i++)
959 x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
960 x3_data[i] += fir_coeffs[n]*x2_data[i];
961 x0_data[i] = x1_data[i];
962 x1_data[i] = x2_data[i];
966 for (
int i = 0; i < num_rows; i++)
968 u_data[i] = x3_data[i];
987 B =
X =
V =
Z = NULL;
993 int _relax_times,
double _relax_weight,
double _omega,
994 int _poly_order,
double _poly_fraction)
1004 B =
X =
V =
Z = NULL;
1013 type =
static_cast<int>(_type);
1039 double a = -1, b, c;
1040 if (!strcmp(name,
"Rectangular")) a = 1.0, b = 0.0, c = 0.0;
1041 if (!strcmp(name,
"Hanning")) a = 0.5, b = 0.5, c = 0.0;
1042 if (!strcmp(name,
"Hamming")) a = 0.54, b = 0.46, c = 0.0;
1043 if (!strcmp(name,
"Blackman")) a = 0.42, b = 0.50, c = 0.08;
1045 mfem_error(
"HypreSmoother::SetWindowByName : name not recognized!");
1061 mfem_error(
"HypreSmoother::SetOperator : not HypreParMatrix!");
1086 A->
Mult(ones, diag);
1099 else if (
type == 1001 ||
type == 1002)
1128 double* window_coeffs =
new double[
poly_order+1];
1129 double* cheby_coeffs =
new double[
poly_order+1];
1137 window_coeffs[i] = a + b*cos(t) +c*cos(2*t);
1141 double theta_pb = acos(1.0 -0.5*k_pb);
1143 cheby_coeffs[0] = (theta_pb +sigma)/M_PI;
1146 double t = i*(theta_pb+sigma);
1147 cheby_coeffs[i] = 2.0*sin(t)/(i*M_PI);
1152 fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
1155 delete[] window_coeffs;
1156 delete[] cheby_coeffs;
1163 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1171 HYPRE_ParCSRDiagScale(NULL, *
A, b, x);
1191 else if (
type == 1002)
1206 hypre_ParCSRRelax(*
A, b,
type,
1211 hypre_ParCSRRelax(*
A, b,
type,
1222 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1228 A -> GetGlobalNumRows(),
1230 A -> GetRowStarts());
1232 A -> GetGlobalNumCols(),
1234 A -> GetColStarts());
1268 :
Solver(_A->Height(), _A->Width())
1279 mfem_error(
"HypreSolver::Mult (...) : HypreParMatrix A is missing");
1297 mfem_error(
"HypreSolver::Mult (...) : HypreParMatrix A is missing");
1303 A -> GetGlobalNumRows(),
1305 A -> GetRowStarts());
1307 A -> GetGlobalNumCols(),
1309 A -> GetColStarts());
1334 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
1336 HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
1341 HYPRE_ParCSRPCGSetTol(pcg_solver, tol);
1346 HYPRE_ParCSRPCGSetMaxIter(pcg_solver, max_iter);
1351 HYPRE_ParCSRPCGSetLogging(pcg_solver, logging);
1356 print_level = print_lvl;
1357 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
1362 HYPRE_ParCSRPCGSetPrecond(pcg_solver,
1370 HYPRE_PCGSetTwoNorm(pcg_solver, 1);
1371 if (res_frequency > 0)
1372 HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
1374 HYPRE_PCGSetResidualTol(pcg_solver, rtol);
1382 double final_res_norm;
1385 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
1389 if (print_level > 0)
1391 time_index = hypre_InitializeTiming(
"PCG Setup");
1392 hypre_BeginTiming(time_index);
1395 HYPRE_ParCSRPCGSetup(pcg_solver, *
A, b, x);
1398 if (print_level > 0)
1400 hypre_EndTiming(time_index);
1401 hypre_PrintTiming(
"Setup phase times", comm);
1402 hypre_FinalizeTiming(time_index);
1403 hypre_ClearTiming();
1407 if (print_level > 0)
1409 time_index = hypre_InitializeTiming(
"PCG Solve");
1410 hypre_BeginTiming(time_index);
1416 HYPRE_ParCSRPCGSolve(pcg_solver, *
A, b, x);
1418 if (print_level > 0)
1420 hypre_EndTiming(time_index);
1421 hypre_PrintTiming(
"Solve phase times", comm);
1422 hypre_FinalizeTiming(time_index);
1423 hypre_ClearTiming();
1425 HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
1426 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
1429 MPI_Comm_rank(comm, &myid);
1433 cout <<
"PCG Iterations = " << num_iterations << endl
1434 <<
"Final PCG Relative Residual Norm = " << final_res_norm
1442 HYPRE_ParCSRPCGDestroy(pcg_solver);
1457 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
1459 HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
1460 HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
1461 HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
1462 HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
1467 HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
1472 HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
1477 HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
1482 HYPRE_ParCSRGMRESSetLogging(gmres_solver, logging);
1487 print_level = print_lvl;
1488 HYPRE_ParCSRGMRESSetPrintLevel(gmres_solver, print_level);
1493 HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
1504 double final_res_norm;
1507 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
1511 if (print_level > 0)
1513 time_index = hypre_InitializeTiming(
"GMRES Setup");
1514 hypre_BeginTiming(time_index);
1517 HYPRE_ParCSRGMRESSetup(gmres_solver, *
A, b, x);
1520 if (print_level > 0)
1522 hypre_EndTiming(time_index);
1523 hypre_PrintTiming(
"Setup phase times", comm);
1524 hypre_FinalizeTiming(time_index);
1525 hypre_ClearTiming();
1529 if (print_level > 0)
1531 time_index = hypre_InitializeTiming(
"GMRES Solve");
1532 hypre_BeginTiming(time_index);
1538 HYPRE_ParCSRGMRESSolve(gmres_solver, *
A, b, x);
1540 if (print_level > 0)
1542 hypre_EndTiming(time_index);
1543 hypre_PrintTiming(
"Solve phase times", comm);
1544 hypre_FinalizeTiming(time_index);
1545 hypre_ClearTiming();
1547 HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
1548 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
1551 MPI_Comm_rank(comm, &myid);
1555 cout <<
"GMRES Iterations = " << num_iterations << endl
1556 <<
"Final GMRES Relative Residual Norm = " << final_res_norm
1564 HYPRE_ParCSRGMRESDestroy(gmres_solver);
1572 int sai_max_levels = 1;
1573 double sai_threshold = 0.1;
1574 double sai_filter = 0.1;
1576 double sai_loadbal = 0.0;
1578 int sai_logging = 1;
1580 HYPRE_ParCSRMatrixGetComm(A, &comm);
1582 HYPRE_ParaSailsCreate(comm, &sai_precond);
1583 HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
1584 HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
1585 HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
1586 HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
1587 HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
1588 HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
1593 HYPRE_ParaSailsSetSym(sai_precond, sym);
1598 HYPRE_ParaSailsDestroy(sai_precond);
1604 int coarsen_type = 10;
1607 int relax_sweeps = 1;
1608 double theta = 0.25;
1609 int interp_type = 6;
1611 int print_level = 1;
1613 HYPRE_BoomerAMGCreate(&amg_precond);
1615 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
1616 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
1617 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
1618 HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
1619 HYPRE_BoomerAMGSetMaxLevels(amg_precond, 25);
1620 HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
1621 HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
1622 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
1623 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
1624 HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
1625 HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
1630 HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
1633 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
1634 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
1639 HYPRE_BoomerAMGDestroy(amg_precond);
1646 int cycle_type = 13;
1649 double rlx_weight = 1.0;
1650 double rlx_omega = 1.0;
1651 int amg_coarsen_type = 10;
1652 int amg_agg_levels = 1;
1653 int amg_rlx_type = 8;
1654 double theta = 0.25;
1655 int amg_interp_type = 6;
1659 if (edge_fespace->
GetNE() > 0)
1663 HYPRE_AMSCreate(&ams);
1665 HYPRE_AMSSetDimension(ams, dim);
1666 HYPRE_AMSSetTol(ams, 0.0);
1667 HYPRE_AMSSetMaxIter(ams, 1);
1668 HYPRE_AMSSetCycleType(ams, cycle_type);
1669 HYPRE_AMSSetPrintLevel(ams, 1);
1683 for (
int i = 0; i < pmesh->
GetNV(); i++)
1685 coord = pmesh -> GetVertex(i);
1686 x_coord(i) = coord[0];
1687 y_coord(i) = coord[1];
1689 z_coord(i) = coord[2];
1696 HYPRE_AMSSetCoordinateVectors(ams, *x, *y, NULL);
1701 HYPRE_AMSSetCoordinateVectors(ams, *x, *y, *z);
1718 HYPRE_AMSSetDiscreteGradient(ams, *G);
1722 Pi = Pix = Piy = Piz = NULL;
1726 if (cycle_type < 10)
1738 if (cycle_type < 10)
1747 Pix = Pi_blocks(0,0);
1748 Piy = Pi_blocks(0,1);
1750 Piz = Pi_blocks(0,2);
1755 HYPRE_AMSSetInterpolations(ams, *Pi, *Pix, *Piy, *Piz);
1757 delete vert_fespace_d;
1760 delete vert_fespace;
1764 HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
1765 HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1766 theta, amg_interp_type, amg_Pmax);
1767 HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1768 theta, amg_interp_type, amg_Pmax);
1773 HYPRE_AMSDestroy(ams);
1789 int cycle_type = 11;
1792 double rlx_weight = 1.0;
1793 double rlx_omega = 1.0;
1794 int amg_coarsen_type = 10;
1795 int amg_agg_levels = 1;
1796 int amg_rlx_type = 6;
1797 double theta = 0.25;
1798 int amg_interp_type = 6;
1800 int ams_cycle_type = 14;
1803 if (face_fespace->
GetNE() > 0)
1806 HYPRE_ADSCreate(&ads);
1808 HYPRE_ADSSetTol(ads, 0.0);
1809 HYPRE_ADSSetMaxIter(ads, 1);
1810 HYPRE_ADSSetCycleType(ads, cycle_type);
1811 HYPRE_ADSSetPrintLevel(ads, 1);
1827 for (
int i = 0; i < pmesh->
GetNV(); i++)
1829 coord = pmesh -> GetVertex(i);
1830 x_coord(i) = coord[0];
1831 y_coord(i) = coord[1];
1832 z_coord(i) = coord[2];
1837 HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
1853 HYPRE_ADSSetDiscreteCurl(ads, *C);
1863 HYPRE_ADSSetDiscreteGradient(ads, *G);
1867 RT_Pi = RT_Pix = RT_Piy = RT_Piz = NULL;
1868 ND_Pi = ND_Pix = ND_Piy = ND_Piz = NULL;
1873 if (ams_cycle_type < 10)
1885 if (ams_cycle_type < 10)
1894 ND_Pix = ND_Pi_blocks(0,0);
1895 ND_Piy = ND_Pi_blocks(0,1);
1896 ND_Piz = ND_Pi_blocks(0,2);
1901 if (cycle_type < 10 && ams_cycle_type > 10)
1903 delete vert_fespace_d;
1907 else if (cycle_type > 10 && ams_cycle_type < 10)
1909 delete vert_fespace_d;
1919 if (cycle_type < 10)
1928 RT_Pix = RT_Pi_blocks(0,0);
1929 RT_Piy = RT_Pi_blocks(0,1);
1930 RT_Piz = RT_Pi_blocks(0,2);
1935 HYPRE_ADSSetInterpolations(ads,
1936 *RT_Pi, *RT_Pix, *RT_Piy, *RT_Piz,
1937 *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz);
1939 delete vert_fespace_d;
1943 delete vert_fespace;
1945 delete edge_fespace;
1948 HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
1949 HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1950 theta, amg_interp_type, amg_Pmax);
1951 HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
1952 amg_rlx_type, theta, amg_interp_type, amg_Pmax);
1957 HYPRE_ADSDestroy(ads);
virtual ~HypreBoomerAMG()
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
int Size() const
Logical size of the array.
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
HypreParVector * X0
FIR Filter Temporary Vectors.
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
int setup_called
Was hypre's Setup function called already?
void SetSize(int s)
Resizes the vector if the new size is different.
HypreParVector * B
Right-hand side and solution vector.
double window_params[3]
Paramters for windowing function of FIR filter.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
HypreADS(HypreParMatrix &A, ParFiniteElementSpace *face_fespace)
void SetWindowByName(const char *window_name)
Convenience function for setting canonical windowing parameters.
HypreParMatrix(hypre_ParCSRMatrix *a)
Converts hypre's format to HypreParMatrix.
int Size() const
Returns the size of the vector.
Abstract parallel finite element space.
bool iterative_mode
If true, use the second argument of Mult as an initial guess.
void Print(const char *fname, int offi=0, int offj=0)
Prints the locally owned rows in parallel.
void SetPrintLevel(int print_lvl)
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
int Size_of_connections() const
double poly_fraction
Fraction of spectrum to smooth for polynomial relaxation.
void SetPrintLevel(int print_lvl)
int poly_scale
Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}.
void AddDomainInterpolator(DiscreteInterpolator *di)
HypreParVector(MPI_Comm comm, int glob_size, int *col)
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's GMRES.
void SetSystemsOptions(int dim)
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
int GetNE() const
Returns number of elements in the mesh.
void SetSymmetry(int sym)
virtual ~HypreParaSails()
void Print(const char *fname)
Prints the locally owned rows in parallel.
int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
MPI_Comm GetComm()
MPI communicator.
double * l1_norms
l1 norms of the rows of A
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
void SetLogging(int logging)
Mesh * GetMesh() const
Returns the mesh.
HyprePCG(HypreParMatrix &_A)
virtual HYPRE_PtrToParSolverFcn SolveFcn() const =0
hypre's internal Solve function
void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0)
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
int Randomize(int seed)
Set random values.
double relax_weight
Damping coefficient (usually <= 1)
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, Array< int > &ess_dof_list, HypreParVector &x, HypreParVector &b)
int Size() const
Returns the number of TYPE I elements.
HypreParMatrix * A
The linear system matrix.
double * GetData() const
Return element data.
HypreAMS(HypreParMatrix &A, ParFiniteElementSpace *edge_fespace)
int * GetI() const
Return the array I.
virtual void SetOperator(const Operator &op)
virtual HYPRE_PtrToParSolverFcn SetupFcn() const =0
hypre's internal Setup function
void SetLogging(int logging)
void SetMaxIter(int max_iter)
int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, hypre_ParVector *f, double max_eig, int poly_order, double *fir_coeffs, hypre_ParVector *u, hypre_ParVector *x0, hypre_ParVector *x1, hypre_ParVector *x2, hypre_ParVector *x3)
void SetSOROptions(double relax_weight, double omega)
Set SOR-related parameters.
Wrapper for hypre's parallel vector class.
int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, hypre_ParVector *f, double lambda, double mu, int N, double max_eig, hypre_ParVector *u, hypre_ParVector *r)
int * GetColStarts() const
Vector * GlobalVector()
Returns the global vector in each processor.
void mfem_error(const char *msg)
int GetOrder(int i) const
Returns the order of the i'th finite element.
int GetGlobalNumCols() const
int GetNumRows() const
Returns the number of rows in the diagonal block of the ParCSRMatrix.
double max_eig_est
Maximal eigenvalue estimate for polynomial smoothing.
void SetDataAndSize(double *d, int s)
HypreParMatrix * Transpose()
Returns the transpose of *this.
void SetFIRCoefficients(double max_eig)
Compute window and Chebyshev coefficients for given polynomial order.
double * fir_coeffs
Combined coefficients for windowing and Chebyshev polynomials.
double InnerProduct(HypreParVector *x, HypreParVector *y)
int * GetRowStarts() const
void SetData(double *_data)
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
void SetTaubinOptions(double lambda, double mu, int iter)
Set parameters for Taubin's lambda-mu method.
HypreParVector * B
Right-hand side and solution vectors.
void SetCommPkg(hypre_ParCSRCommPkg *comm_pkg)
int GetNumCols() const
Returns the number of columns in the diagonal block of the ParCSRMatrix.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
void SetPolyOptions(int poly_order, double poly_fraction)
Set parameters for polynomial smoothing.
int relax_times
Number of relaxation sweeps.
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
int * GetTrueDofOffsets()
Abstract class for hypre's solvers and preconditioners.
int GetGlobalNumRows() const
HypreParVector & operator=(double d)
Set constant values.
Arbitrary order H(curl)-conforming Nedelec finite elements.
void GetDiag(Vector &diag)
Get the diagonal of the matrix.
HypreGMRES(HypreParMatrix &_A)
virtual ~HypreParMatrix()
Calls hypre's destroy function.
void GetParBlocks(Array2D< HypreParMatrix * > &blocks) const
hypre_ParCSRMatrix * StealData()
Changes the ownership of the the matrix.
HypreParMatrix * A
The linear system matrix.
Arbitrary order H1-conforming (continuous) finite elements.
void SetMaxIter(int max_iter)
double omega
SOR parameter (usually in (0,2))
Class for parallel grid function.
HypreParMatrix * ParallelAssemble(SparseMatrix *m)
Wrapper for hypre's ParCSR matrix class.
~HypreParVector()
Calls hypre's destroy function.
virtual void Assemble(int skip_zeros=1)
int * GetJ() const
Return the array J.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Class for parallel meshes.
void Read(MPI_Comm comm, const char *fname)
Reads the matrix from a file.
void ParallelAverage(Vector &tv) const
Returns the vector averaged on the true dofs.
HypreParVector * V
Temporary vectors.
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
HypreParaSails(HypreParMatrix &A)
int poly_order
Order of the smoothing polynomial.
double lambda
Taubin's lambda-mu method parameters.
HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B)
Returns the matrix A * B.
HypreBoomerAMG(HypreParMatrix &A)