本文整理汇总了C++中MatGetRow函数的典型用法代码示例。如果您正苦于以下问题:C++ MatGetRow函数的具体用法?C++ MatGetRow怎么用?C++ MatGetRow使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MatGetRow函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: MatAXPY_Basic
PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
{
PetscInt i,start,end,j,ncols,m,n;
PetscErrorCode ierr;
const PetscInt *row;
PetscScalar *val;
const PetscScalar *vals;
PetscFunctionBegin;
ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr);
if (a == 1.0) {
for (i = start; i < end; i++) {
ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr);
ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
}
} else {
ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&val);CHKERRQ(ierr);
for (i=start; i<end; i++) {
ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
for (j=0; j<ncols; j++) {
val[j] = a*vals[j];
}
ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr);
ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
}
ierr = PetscFree(val);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:feelpp,项目名称:debian-petsc,代码行数:33,代码来源:axpy.c
示例2: testBSplineSetD2R1Mat
int testBSplineSetD2R1Mat() {
PrintTimeStamp(PETSC_COMM_SELF, "D2", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
int order = 3;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
BSSSetUp(bss);
// compute matrix
Mat M; BSSCreateR1Mat(bss, &M); BSSD2R1Mat(bss, M);
// value check
const PetscScalar *row;
MatGetRow(M, 1, NULL, NULL, &row);
ASSERT_DOUBLE_NEAR(0.1666666666, row[0], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(-1.0, row[1], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.3333333333, row[2], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.16666666666, row[3], pow(10.0, -10.0));
MatRestoreRow(M, 1, NULL, NULL, &row);
MatGetRow(M, 2, NULL, NULL, &row);
ASSERT_DOUBLE_NEAR(0.1666666666666, row[0], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.3333333333333, row[1], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(-1.0, row[2] , pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.3333333333333, row[3], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.1666666666666, row[4], pow(10.0, -10.0));
MatRestoreRow(M, 2, NULL, NULL, &row);
// Finalize
MatDestroy(&M);
BSSDestroy(&bss);
return 0;
}
开发者ID:ReiMatsuzaki,项目名称:rescol,代码行数:35,代码来源:test_bspline.c
示例3: TEST
TEST(schur, complement) {
Profiler::initialize();
IS set;
// vytvorit rozdeleni bloku na procesory ve tvaru "part" (tj. indexy prvnich radku na procesorech)
int np, rank;
MPI_Comm_size(PETSC_COMM_WORLD, &np);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
Distribution ds(block_size, MPI_COMM_WORLD);
Distribution block_ds(block_count, MPI_COMM_WORLD);
Distribution all_ds(block_size + block_count, MPI_COMM_WORLD);
/*if (rank == 0) {
cout << all_ds;
cout << ds;
cout << block_ds;
}*/
ISCreateStride(PETSC_COMM_WORLD, ds.lsize(), all_ds.begin(), 1, &set);
ISView(set, PETSC_VIEWER_STDOUT_WORLD);
// volat s lokalni velkosti = pocet radku na lokalnim proc.
SchurComplementTest * schurComplement = new SchurComplementTest(set, &all_ds);
schurComplement->set_solution(NULL);
schurComplement->set_positive_definite();
schurComplement->start_allocation();
schurComplement->fill_matrix( rank, ds, block_ds); // preallocate matrix
schurComplement->start_add_assembly();
schurComplement->fill_matrix( rank, ds, block_ds); // fill matrix
schurComplement->finish_assembly();
MatView(*(schurComplement->get_matrix()),PETSC_VIEWER_STDOUT_WORLD);
LinSys * lin_sys = new LinSysPetscTest( schurComplement->make_complement_distribution() );
schurComplement->set_complement( (LinSys_PETSC *)lin_sys );
schurComplement->solve();
// test of computed values
{
PetscInt ncols;
const PetscInt *cols;
const PetscScalar *vals;
for (unsigned int i=0; i<block_size; i++) {
MatGetRow(schurComplement->get_a_inv(), i + rank*block_size, &ncols, &cols, &vals);
EXPECT_FLOAT_EQ( (1.0 / (double)(rank + 2)), vals[i] );
MatRestoreRow(schurComplement->get_a_inv(), i + rank*block_size, &ncols, &cols, &vals);
}
MatGetRow(*(schurComplement->get_system()->get_matrix()), rank, &ncols, &cols, &vals);
EXPECT_FLOAT_EQ( ((double)block_size / (double)(rank + 2)), vals[0] );
MatRestoreRow(*(schurComplement->get_system()->get_matrix()), rank, &ncols, &cols, &vals);
}
}
开发者ID:jbrezmorf,项目名称:flow123d,代码行数:53,代码来源:schur_compl_test.cpp
示例4: testBSplineSetSR1Mat
int testBSplineSetSR1Mat() {
PrintTimeStamp(PETSC_COMM_SELF, "S", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
int order = 3;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
BSSSetUp(bss);
// compute S matrix
Mat S;
PetscErrorCode ierr;
ierr = BSSCreateR1Mat(bss, &S);
ierr = BSSSR1Mat(bss, S); CHKERRQ(ierr);
// Get structure and check
int m, n;
ierr = MatGetSize(S, &m, &n); CHKERRQ(ierr);
ASSERT_EQ(5, m);
ASSERT_EQ(5, n);
// value check
const PetscScalar *row;
PetscInt ncols;
const PetscInt *cols;
ierr = MatGetRow(S, 0, &ncols, &cols, &row); CHKERRQ(ierr);
ASSERT_EQ(3, ncols);
ASSERT_EQ(0, cols[0]);
ASSERT_DOUBLE_EQ(1.0/3.0, row[0]);
ASSERT_EQ(1, cols[1]);
ASSERT_DOUBLE_NEAR(0.2083333333, row[1], pow(10.0, -10.0));
ASSERT_EQ(2, cols[2]);
ASSERT_DOUBLE_NEAR(0.0083333333, row[2], pow(10.0, -10.0));
ierr = MatRestoreRow(S, 0, &ncols, &cols, &row); CHKERRQ(ierr);
ierr = MatGetRow(S, 2, &ncols, &cols, &row); CHKERRQ(ierr);
ASSERT_EQ(5, ncols);
for(int i = 0; i < 5; i++)
ASSERT_EQ(i, cols[i]);
ASSERT_DOUBLE_NEAR(0.0083333333333, row[0], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.2166666666666, row[1], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.55, row[2], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.2166666666666, row[3], pow(10.0, -10.0));
ASSERT_DOUBLE_NEAR(0.0083333333, row[4], pow(10.0, -10.0));
ierr = MatRestoreRow(S, 3, &ncols, &cols, &row); CHKERRQ(ierr);
// Finalize
MatDestroy(&S);
BSSDestroy(&bss);
return 0;
}
开发者ID:ReiMatsuzaki,项目名称:rescol,代码行数:52,代码来源:test_bspline.c
示例5: MatCreate
/*@
MatChop - Set all values in the matrix less than the tolerance to zero
Input Parameters:
+ A - The matrix
- tol - The zero tolerance
Output Parameters:
. A - The chopped matrix
Level: intermediate
.seealso: MatCreate(), MatZeroEntries()
@*/
PetscErrorCode MatChop(Mat A, PetscReal tol)
{
PetscScalar *newVals;
PetscInt *newCols;
PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = MatGetOwnershipRange(A, &rStart, &rEnd);
CHKERRQ(ierr);
for (r = rStart; r < rEnd; ++r) {
PetscInt ncols;
ierr = MatGetRow(A, r, &ncols, NULL, NULL);
CHKERRQ(ierr);
colMax = PetscMax(colMax, ncols);
CHKERRQ(ierr);
ierr = MatRestoreRow(A, r, &ncols, NULL, NULL);
CHKERRQ(ierr);
}
numRows = rEnd - rStart;
ierr = MPI_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));
CHKERRQ(ierr);
ierr = PetscMalloc2(colMax,&newCols,colMax,&newVals);
CHKERRQ(ierr);
for (r = rStart; r < rStart+maxRows; ++r) {
const PetscScalar *vals;
const PetscInt *cols;
PetscInt ncols, newcols, c;
if (r < rEnd) {
ierr = MatGetRow(A, r, &ncols, &cols, &vals);
CHKERRQ(ierr);
for (c = 0; c < ncols; ++c) {
newCols[c] = cols[c];
newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
}
newcols = ncols;
ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);
CHKERRQ(ierr);
ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
}
ierr = PetscFree2(newCols,newVals);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:pombredanne,项目名称:petsc,代码行数:66,代码来源:axpy.c
示例6: MatGetRow
int Epetra_PETScAIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const
{
int globalRow = PetscRowStart_ + Row;
MatGetRow(Amat_, globalRow, &NumEntries, PETSC_NULL, PETSC_NULL);
MatRestoreRow(Amat_,globalRow,&NumEntries, PETSC_NULL, PETSC_NULL);
return(0);
}
开发者ID:00liujj,项目名称:trilinos,代码行数:7,代码来源:EpetraExt_PETScAIJMatrix.cpp
示例7: MatGetRow
void PETSc::Get_row_of_A(PetscInt i, PetscInt *ncol, PetscInt **cols, double **row)
{
ierr = MatGetRow(A,i,ncol,(const PetscInt**)cols,
(const PetscScalar**)row);
ierr = MatRestoreRow(A,i,ncol,(const PetscInt**)cols,
(const PetscScalar**)row);
}
开发者ID:tonyguo1,项目名称:LSGFD,代码行数:7,代码来源:solver_petsc.cpp
示例8: MatGetRow
void PetscSparseStorage::fromGetSpRow( int row, int col,
double A[], int /*lenA */,
int jcolA[], int& nnz,
int colExtent, int& info )
{
const double * B;
const int * jcolB;
int k, ierr;
int m, n;
this->getSize(m,n);
ierr = MatGetRow(M, row, &nnz, &jcolB, &B ); assert( ierr == 0 );
if( col == 0 && colExtent == n ) {
// We want the whole row, just do a simple copy
for( k = 0; k < nnz; k++ ) {
A[k] = B[k];
jcolA[k] = jcolB[k];
}
} else { // Copy only those in range
int i = 0;
for( k = 0; k < nnz; k++ ) {
if( jcolB[k] >= col && jcolB[k] < col + colExtent ) {
A[i] = B[k];
jcolA[i] = jcolB[k];
i++;
}
}
}
ierr = MatRestoreRow( M, row, &nnz, &jcolB, &B ); assert( ierr == 0 );
info = 0;
}
开发者ID:aig-lchion,项目名称:OOQP,代码行数:31,代码来源:PetscSparseStorage.C
示例9: testBSplineSetENMatR1Mat
int testBSplineSetENMatR1Mat() {
PrintTimeStamp(PETSC_COMM_SELF, "EN", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
int order = 3;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
BSSSetUp(bss);
// compute matrix
Mat M; BSSCreateR1Mat(bss, &M); BSSENR1Mat(bss, 2, 0.7, M);
// value check
const PetscScalar *row;
MatGetRow(M, 4, NULL, NULL, &row);
ASSERT_DOUBLE_NEAR(0.0000964408077, row[0], pow(10.0, -8));
ASSERT_DOUBLE_NEAR(0.00168615366, row[1], pow(10.0, -8));
ASSERT_DOUBLE_NEAR(0.00211136230, row[2], pow(10.0, -8));
MatRestoreRow(M, 4, NULL, NULL, &row);
// Finalize
MatDestroy(&M);
BSSDestroy(&bss);
return 0;
}
开发者ID:ReiMatsuzaki,项目名称:rescol,代码行数:26,代码来源:test_bspline.c
示例10: MatDumpSPAI
/*
MatDumpSPAI - Dumps a PETSc matrix to a file in an ASCII format
suitable for the SPAI code of Stephen Barnard to solve. This routine
is simply here to allow testing of matrices directly with the SPAI
code, rather then through the PETSc interface.
*/
PetscErrorCode MatDumpSPAI(Mat A,FILE *file)
{
const PetscScalar *vals;
PetscErrorCode ierr;
int i,j,n,size,nz;
const int *cols;
MPI_Comm comm;
PetscObjectGetComm((PetscObject)A,&comm);
MPI_Comm_size(comm,&size);
if (size > 1) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Only single processor dumps");
ierr = MatGetSize(A,&n,&n);CHKERRQ(ierr);
/* print the matrix */
fprintf(file,"%d\n",n);
for (i=0; i<n; i++) {
ierr = MatGetRow(A,i,&nz,&cols,&vals);CHKERRQ(ierr);
for (j=0; j<nz; j++) {
fprintf(file,"%d %d %16.14e\n",i+1,cols[j]+1,vals[j]);
}
ierr = MatRestoreRow(A,i,&nz,&cols,&vals);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:34,代码来源:dspai.c
示例11: DMPlexCreate
/*@
MatComputeBandwidth - Calculate the full bandwidth of the matrix, meaning the width 2k+1 where k diagonals on either side are sufficient to contain all the matrix nonzeros.
Collective on Mat
Input Parameters:
+ A - The Mat
- fraction - An optional percentage of the Frobenius norm of the matrix that the bandwidth should enclose
Output Parameter:
. bw - The matrix bandwidth
Level: beginner
.seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
@*/
PetscErrorCode MatComputeBandwidth(Mat A, PetscReal fraction, PetscInt *bw)
{
PetscInt lbw[2] = {0, 0}, gbw[2];
PetscInt rStart, rEnd, r;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
PetscValidLogicalCollectiveReal(A,fraction,2);
PetscValidPointer(bw, 3);
if ((fraction > 0.0) && (fraction < 1.0)) SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "We do not yet support a fractional bandwidth");
ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr);
for (r = rStart; r < rEnd; ++r) {
const PetscInt *cols;
PetscInt ncols;
ierr = MatGetRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
if (ncols) {
lbw[0] = PetscMax(lbw[0], r - cols[0]);
lbw[1] = PetscMax(lbw[1], cols[ncols-1] - r);
}
ierr = MatRestoreRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
}
ierr = MPI_Allreduce(lbw, gbw, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) A));CHKERRQ(ierr);
*bw = 2*PetscMax(gbw[0], gbw[1]) + 1;
PetscFunctionReturn(0);
}
开发者ID:pombredanne,项目名称:petsc,代码行数:43,代码来源:bandwidth.c
示例12: MatGetDiagonalHermitian_Normal
PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
{
Mat_Normal *Na = (Mat_Normal*)N->data;
Mat A = Na->A;
PetscErrorCode ierr;
PetscInt i,j,rstart,rend,nnz;
const PetscInt *cols;
PetscScalar *diag,*work,*values;
const PetscScalar *mvalues;
PetscFunctionBegin;
ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr);
ierr = PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
for (i=rstart; i<rend; i++) {
ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
for (j=0; j<nnz; j++) {
work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
}
ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
}
ierr = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr);
rstart = N->cmap->rstart;
rend = N->cmap->rend;
ierr = VecGetArray(v,&values);CHKERRQ(ierr);
ierr = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr);
ierr = VecRestoreArray(v,&values);CHKERRQ(ierr);
ierr = PetscFree2(diag,work);CHKERRQ(ierr);
ierr = VecScale(v,Na->scale);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:masa-ito,项目名称:PETScToPoisson,代码行数:31,代码来源:normmh.c
示例13: times
void
PetscSparseMtrx :: times(const FloatMatrix &B, FloatMatrix &answer) const
{
if ( this->giveNumberOfColumns() != B.giveNumberOfRows() ) {
OOFEM_ERROR("Dimension mismatch");
}
#ifdef __PARALLEL_MODE
if ( emodel->isParallel() ) {
OOFEM_ERROR("PetscSparseMtrx :: times - Not implemented");
}
#endif
// I'm opting to work with a set of vectors, as i think it might be faster and more robust. / Mikael
int nr = this->giveNumberOfRows();
int nc = B.giveNumberOfColumns();
answer.resize(nr, nc);
double *aptr = answer.givePointer();
#if 0
// Approach using several vectors. Not sure if it is optimal, but it includes petsc calls which i suspect are inefficient. / Mikael
// UNTESTED!
Vec globX, globY;
VecCreate(PETSC_COMM_SELF, &globY);
VecSetType(globY, VECSEQ);
VecSetSizes(globY, PETSC_DECIDE, nr);
int nrB = B.giveNumberOfRows();
for (int k = 0; k < nc; k++) {
double colVals[nrB];
for (int i = 0; i < nrB; i++) colVals[i] = B(i,k); // B.copyColumn(Bk,k);
VecCreateSeqWithArray(PETSC_COMM_SELF, nrB, colVals, &globX);
MatMult(this->mtrx, globX, globY );
double *ptr;
VecGetArray(globY, &ptr);
for (int i = 0; i < nr; i++) *aptr++ = ptr[i]; // answer.setColumn(Ak,k);
VecRestoreArray(globY, &ptr);
VecDestroy(globX);
}
VecDestroy(globY);
#endif
Mat globB, globC;
MatCreateSeqDense(PETSC_COMM_SELF, B.giveNumberOfRows(), B.giveNumberOfColumns(), B.givePointer(), & globB);
MatMatMult(this->mtrx, globB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, & globC);
const double *vals;
for ( int r = 0; r < nr; r++ ) {
MatGetRow(globC, r, NULL, NULL, & vals);
for ( int i = 0, i2 = r; i < nc; i++, i2 += nr ) {
aptr [ i2 ] = vals [ i ];
}
MatRestoreRow(globC, r, NULL, NULL, & vals);
}
MatDestroy(&globB);
MatDestroy(&globC);
}
开发者ID:MartinFagerstrom,项目名称:oofem,代码行数:56,代码来源:petscsparsemtrx.C
示例14: MatGetRow_SMF
PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
{
PetscErrorCode ierr;
MatSubMatFreeCtx ctx;
PetscFunctionBegin;
ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
ierr = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:tom-klotz,项目名称:petsc,代码行数:10,代码来源:submatfree.c
示例15: PCGAMGOptProl_Classical_Jacobi
PetscErrorCode PCGAMGOptProl_Classical_Jacobi(PC pc,const Mat A,Mat *P)
{
PetscErrorCode ierr;
PetscInt f,s,n,cf,cs,i,idx;
PetscInt *coarserows;
PetscInt ncols;
const PetscInt *pcols;
const PetscScalar *pvals;
Mat Pnew;
Vec diag;
PC_MG *mg = (PC_MG*)pc->data;
PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
PetscFunctionBegin;
if (cls->nsmooths == 0) {
ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
ierr = MatGetOwnershipRange(*P,&s,&f);CHKERRQ(ierr);
n = f-s;
ierr = MatGetOwnershipRangeColumn(*P,&cs,&cf);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscInt)*n,&coarserows);CHKERRQ(ierr);
/* identify the rows corresponding to coarse unknowns */
idx = 0;
for (i=s;i<f;i++) {
ierr = MatGetRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
/* assume, for now, that it's a coarse unknown if it has a single unit entry */
if (ncols == 1) {
if (pvals[0] == 1.) {
coarserows[idx] = i;
idx++;
}
}
ierr = MatRestoreRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
}
ierr = MatGetVecs(A,&diag,0);CHKERRQ(ierr);
ierr = MatGetDiagonal(A,diag);CHKERRQ(ierr);
ierr = VecReciprocal(diag);CHKERRQ(ierr);
for (i=0;i<cls->nsmooths;i++) {
ierr = MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);CHKERRQ(ierr);
ierr = MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);CHKERRQ(ierr);
ierr = MatDiagonalScale(Pnew,diag,0);CHKERRQ(ierr);
ierr = MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = MatDestroy(P);CHKERRQ(ierr);
*P = Pnew;
Pnew = NULL;
}
ierr = VecDestroy(&diag);CHKERRQ(ierr);
ierr = PetscFree(coarserows);CHKERRQ(ierr);
ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:ZJLi2013,项目名称:petsc,代码行数:54,代码来源:classical.c
示例16: MatConvert_MPIAIJ_MPISBAIJ
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "MatConvert_MPIAIJ_MPISBAIJ"
PetscErrorCode MatConvert_MPIAIJ_MPISBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
{
PetscErrorCode ierr;
Mat M;
Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)A->data;
Mat_SeqAIJ *Aa = (Mat_SeqAIJ*)mpimat->A->data,*Ba = (Mat_SeqAIJ*)mpimat->B->data;
PetscInt *d_nnz,*o_nnz;
PetscInt i,j,nz;
PetscInt m,n,lm,ln;
PetscInt rstart,rend;
const PetscScalar *vwork;
const PetscInt *cwork;
PetscFunctionBegin;
if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
ierr = MatGetLocalSize(A,&lm,&ln);CHKERRQ(ierr);
ierr = PetscMalloc2(lm,PetscInt,&d_nnz,lm,PetscInt,&o_nnz);CHKERRQ(ierr);
ierr = MatMarkDiagonal_SeqAIJ(mpimat->A);CHKERRQ(ierr);
for (i=0;i<lm;i++){
d_nnz[i] = Aa->i[i+1] - Aa->diag[i];
o_nnz[i] = Ba->i[i+1] - Ba->i[i];
}
ierr = MatCreate(((PetscObject)A)->comm,&M);CHKERRQ(ierr);
ierr = MatSetSizes(M,lm,ln,m,n);CHKERRQ(ierr);
ierr = MatSetType(M,MATMPISBAIJ);CHKERRQ(ierr);
ierr = MatSeqSBAIJSetPreallocation(M,1,0,d_nnz);CHKERRQ(ierr);
ierr = MatMPISBAIJSetPreallocation(M,1,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
for (i=rstart;i<rend;i++){
ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
j = 0;
while (cwork[j] < i){ j++; nz--;}
ierr = MatSetValues(M,1,&i,nz,cwork+j,vwork+j,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (reuse == MAT_REUSE_MATRIX) {
ierr = MatHeaderReplace(A,M);CHKERRQ(ierr);
} else {
*newmat = M;
}
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:54,代码来源:mpiaijsbaij.c
示例17: build_strength_matrix
void
build_strength_matrix(Mat A, PetscReal theta, Mat* strength) {
//Variables for the new matrix structure
std::vector<PetscInt> rows, cols;
int cursor = 0;
PetscInt start;
PetscInt end;
MatGetOwnershipRange((A), &start, &end);
for (int row=start; row<end; row++) {
PetscInt ncols;
const PetscInt *col_indx;
const PetscScalar *col_value;
MatGetRow((A), row, &ncols, &col_indx, &col_value);
rows.push_back(cursor);
// First, find the threshhold for this row
PetscScalar strong_threshhold = -col_value[0];
for (int ii=0; ii<ncols; ii++) {
if (-col_value[ii] > strong_threshhold) {
strong_threshhold = -col_value[ii];
}
}
strong_threshhold *= theta;
//if the threshold is negative, assume that this row only has a diagonal entry and skip the row
if (strong_threshhold > 0) {
for (int ii=0; ii<ncols; ii++) {
if (-col_value[ii] >= strong_threshhold) {
cols.push_back(col_indx[ii]);
cursor++;
}
}
}
MatRestoreRow((A), row, &ncols, &col_indx, &col_value);
}
rows.push_back(cursor);
MatGetOwnershipRange((A), &start, &end);
std::vector<PetscScalar> data(cols.size());
//TODO: control for cases where cols and rows are split differently
//TODO: replace this PETSC_COMM_WORLD so the strength matrix is using the same communicator as the original matrix.
MatCreate(PETSC_COMM_WORLD, strength);
MatSetSizes(*strength, end-start, end-start, PETSC_DETERMINE, PETSC_DETERMINE);
MatSetType(*strength,MATMPIAIJ);
MatMPIAIJSetPreallocationCSR(*strength,&rows[0],&cols[0],&data[0]);
//TODO: the above code is a work around for a bug in the following function call:
//MatCreateMPIAIJWithArrays(PETSC_COMM_WORLD, end-start, end-start, PETSC_DETERMINE, PETSC_DETERMINE, &rows[0], &cols[0], &data[0], strength);
}
开发者ID:rblake,项目名称:petsc_amg,代码行数:52,代码来源:mglib.c
示例18: MatConvert_MPIBAIJ_MPISBAIJ
PETSC_EXTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
{
PetscErrorCode ierr;
Mat M;
Mat_MPIBAIJ *mpimat = (Mat_MPIBAIJ*)A->data;
Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ*)mpimat->A->data,*Ba = (Mat_SeqBAIJ*)mpimat->B->data;
PetscInt *d_nnz,*o_nnz;
PetscInt i,j,nz;
PetscInt m,n,lm,ln;
PetscInt rstart,rend;
const PetscScalar *vwork;
const PetscInt *cwork;
PetscInt bs = A->rmap->bs;
PetscFunctionBegin;
ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
ierr = MatGetLocalSize(A,&lm,&ln);CHKERRQ(ierr);
ierr = PetscMalloc2(lm/bs,&d_nnz,lm/bs,&o_nnz);CHKERRQ(ierr);
ierr = MatMarkDiagonal_SeqBAIJ(mpimat->A);CHKERRQ(ierr);
for (i=0; i<lm/bs; i++) {
d_nnz[i] = Aa->i[i+1] - Aa->diag[i];
o_nnz[i] = Ba->i[i+1] - Ba->i[i];
}
ierr = MatCreate(PetscObjectComm((PetscObject)A),&M);CHKERRQ(ierr);
ierr = MatSetSizes(M,lm,ln,m,n);CHKERRQ(ierr);
ierr = MatSetType(M,MATMPISBAIJ);CHKERRQ(ierr);
ierr = MatSeqSBAIJSetPreallocation(M,bs,0,d_nnz);CHKERRQ(ierr);
ierr = MatMPISBAIJSetPreallocation(M,bs,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
ierr = MatSetOption(M,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);
for (i=rstart; i<rend; i++) {
ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
j = 0;
ierr = MatSetValues(M,1,&i,nz,cwork+j,vwork+j,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (reuse == MAT_REUSE_MATRIX) {
ierr = MatHeaderReplace(A,&M);CHKERRQ(ierr);
} else {
*newmat = M;
}
PetscFunctionReturn(0);
}
开发者ID:pombredanne,项目名称:petsc,代码行数:51,代码来源:mpiaijsbaij.c
示例19: EuclidGetRow
void EuclidGetRow(void *Ain, HYPRE_Int row, HYPRE_Int *len, HYPRE_Int **ind, double **val)
{
START_FUNC_DH
Mat A = Ain;
HYPRE_Int ierr;
ierr = MatGetRow(A, row, len, ind, val);
if (ierr) {
hypre_sprintf(msgBuf_dh, "PETSc's MatGetRow bombed for row= %i", row);
SET_V_ERROR(msgBuf_dh);
}
END_FUNC_DH
}
开发者ID:tpatki,项目名称:rapl-old-data,代码行数:14,代码来源:getRow_dh.c
示例20: PreallocateJacobian
PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
{
PetscErrorCode ierr;
PetscInt *d_nnz;
PetscInt i,idx=0,start=0;
PetscFunctionBegin;
ierr = PetscMalloc1(user->neqs_pgrid,&d_nnz);CHKERRQ(ierr);
for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0;
/* Generator subsystem */
for (i=0; i < ngen; i++) {
d_nnz[idx] += 3;
d_nnz[idx+1] += 2;
d_nnz[idx+2] += 2;
d_nnz[idx+3] += 5;
d_nnz[idx+4] += 6;
d_nnz[idx+5] += 6;
d_nnz[user->neqs_gen+2*gbus[i]] += 3;
d_nnz[user->neqs_gen+2*gbus[i]+1] += 3;
d_nnz[idx+6] += 2;
d_nnz[idx+7] += 2;
d_nnz[idx+8] += 5;
idx = idx + 9;
}
start = user->neqs_gen;
PetscInt ncols;
for (i=0; i < nbus; i++) {
ierr = MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr);
d_nnz[start+2*i] += ncols;
d_nnz[start+2*i+1] += ncols;
ierr = MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);CHKERRQ(ierr);
}
ierr = MatSeqAIJSetPreallocation(J,0,d_nnz);CHKERRQ(ierr);
ierr = PetscFree(d_nnz);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:00liujj,项目名称:petsc,代码行数:45,代码来源:ex9bus.c
注:本文中的MatGetRow函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论