[petsc-users] Outer Product of two vectors
Barry Smith
bsmith at mcs.anl.gov
Wed Oct 6 14:42:50 CDT 2010
I think VecScatterCreateToAll() is what you need. Something like
VecScatterCreateToAll(x,&ctx,&xeveryone);
VecScatterBegin(ctx,x,everyone,INSERT_VALUES);
VecGetArray(xeveryone,&xvalues);
Now every process has all the values of x in the array xvalues
On Oct 6, 2010, at 2:38 PM, Harsh Sharma wrote:
> Hi,
>
> I noticed that PETSc doesn't have a routine to compute
> the outer product of two Vec objects.
>
> Specifically, I want to compute the outer product for two
> Vec objects that are large enough to require storage on
> multiple processors.
>
> If I write my own routine in such a way that it performs
> local computations of the entries of the outer product
> matrix (called OP from hereon), say using
> VecGetOwnershipRange to retrieve local entries of the
> two vectors (Vec a and Vec b) and then computing local
> blocks of OP, then there will be parts of OP that will not
> be computed.
>
> Say we have 3 processors, and
> Vec a is (2,4,8,16,32), Vec b is (1, 1/2, 1/4). Then the
> above approach will produce OP as
>
> 2 - -
> 4 - -
> - 4 -
> - 8 -
> - - 8
>
> instead of
>
> 2 1 0.5
> 4 2 1
> 8 4 2
> 16 8 4
> 32 16 8
>
> So, I chose to implement my routine in this way instead >>
> I would convert the "row" vector {vector b in a*transpose(b) }
> into an array of PetscScalar values, using PetscMalloc.
> Then, I run a for-loop traversing this array-version of vector b
> where in each iteration of the loop, the processor-specific
> part of vector a is scaled with the element of b that is being
> accessed in the loop. Then MatSetValues is called to set
> the scaled part-of-vector-a in the right locations in OP.
>
> However, this approach is also producing the same result as
> above. My guess is that the array-version of Vec b created
> using PetscMalloc is not a globally visible array that is being accessed
> by each processor -- instead, each processor is creating its own
> version of array-vector-b and then again computing only parts of
> OP. Further, I think this probably has to do with the fact that
> the routine PetscMalloc is not collective.
>
> How do I get an array of PetscScalar values that is not
> processor-specific but is visible to all processors? If this
> cannot be done, how do I go about computing the outer
> product of two vectors?
>
> I am appending my code (routine to compute outer product and
> the "main" function) and sample output with this mail.
>
> Thanks very much,
> Harsh
>
> ---- code and sample output ----
>
> static char helpMsg[] = "\nComputes outer product of two vectors.\n";
> #include "petscmat.h"
>
> // function to compute outer-product matrix of two vectors
> PetscErrorCode MyMPIVecOuterProd(Mat OP, Vec a, Vec b, InsertMode addV)
> {
> /*
> for vectors a and b, computes a.transpose(b) and
> adds/stores the resulting outer-product matrix
> to/in the matrix OP
> */
>
> PetscInt nRows = 0; /* number of rows of OP matrix */
> PetscInt nCols = 0; /* number of columns of OP matrix */
> PetscInt nA = 0; /* length of vector a */
> PetscInt nB = 0; /* length of vector b */
>
> PetscScalar *locAVals; /* array to hold local vector a values */
> PetscScalar *locSAVals; /* array to hold scaled local vector a values */
> PetscScalar *locBVals; /* array to hold local vector b values */
> PetscScalar *bArr; /* array to hold the entire vector b */
>
> PetscInt aLow,aHigh,bLow,bHigh; /* local index-range limits */
> PetscInt ia,ib; /* for-loop index variables for a and b vectors */
> PetscInt * locRowIdxOP; /* locally-set OP column indices for MatSetValues */
>
> PetscScalar curBVal; /* value of vector b's component with which to scale vector a */
>
> /* get the dimension of vector a */
> VecGetSize(a,&nA);
> /* get the dimension of vector b */
> VecGetSize(b,&nB);
>
> /* get the dimensions of outer-product matrix */
> MatGetSize(OP,&nRows,&nCols);
>
> /* check for dimensional compatibility */
> if ((nRows != nA) || (nCols != nB)) {
> SETERRQ(1,"Error: MyMPIVecOuterProd: Dimensional Incompatibility!");
> return(1);
> }
>
> /* --------------------------------------------- */
> /* first, convert vector b into array of scalars */
> /* --------------------------------------------- */
>
> /* allocate memory for array-representation of vector b */
> PetscMalloc((nB)*sizeof(PetscScalar),&bArr);
>
> /* do local assignment from vector b values to array-representation */
> /* first, obtain local range of vector b */
> VecGetOwnershipRange(b,&bLow,&bHigh); // bHigh is one more than highest local index
> /* then, obtain pointer to local elements of vector b */
> VecGetArray(b,&locBVals);
> /* then, assign local values of vector b to corresponding locations in bArr */
> for (ib = bLow; ib < bHigh; ib++) {
> *(bArr + ib) = *(locBVals + ib - bLow);
> } // end of b for loop
> /* finally, restore local elements of vector b */
> VecRestoreArray(b,&locBVals);
>
> /* ------------------------------------------------- */
> /* next, scale local values of vector a and add them */
> /* to corresponding locations in the OP matrix */
> /* ------------------------------------------------- */
>
> /* first, obtain local range of vector a */
> VecGetOwnershipRange(a,&aLow,&aHigh); // aHigh is one more than highest local index
> /* then, obtain pointer to local elements of vector a */
> VecGetArray(a,&locAVals);
> /* allocate memory for local array-of-scaled-vector-a-values */
> PetscMalloc((aHigh-aLow)*sizeof(PetscScalar),&locSAVals);
> /* allocate memory for locally-set OP row indices */
> PetscMalloc((aHigh-aLow)*sizeof(PetscInt),&locRowIdxOP);
>
> /* set locally-set OP row indices */
> for (ia = 0; ia < aHigh-aLow; ia++) {
> *(locRowIdxOP + ia) = ia + aLow;
> } // end of for : set locally-set OP row indices
>
> /* next, for each component of vector b (bArr), scale local vector a values */
> /* and set them up in the corresponding locations in the OP matrix */
> for (ib = 0; ib < nB; ib++) {
> /* get component of vector b to scale with */
> curBVal = *(bArr + ib);
> /* scale vector a local values */
> for (ia = 0; ia < aHigh-aLow; ia++) {
> *(locSAVals + ia) = (*(locAVals + ia)) * curBVal;
> } // end of for: scale local vector a values
> /* set scaled values in appropriate locations in OP */
> MatSetValues(OP,aHigh-aLow,locRowIdxOP,1,&ib,locSAVals,addV)
> } // end of for: set scaled values in OP matrix
>
> /* next, restore local elements of vector a */
> VecRestoreArray(a,&locAVals);
> /* free memory for local row indices of OP */
> PetscFree(locRowIdxOP);
> /* free memory for local scaled vector a values */
> PetscFree(locSAVals);
> /* free memory for array representation of vector b */
> PetscFree(bArr);
>
> /* ------------------------------- */
> /* finally, assemble the OP matrix */
> /* ------------------------------- */
> MatAssemblyBegin(OP,MAT_FINAL_ASSEMBLY);
> MatAssemblyEnd(OP,MAT_FINAL_ASSEMBLY);
>
> return(0);
> }
>
> int main (int argc, char **argv) {
> PetscInt n1=5; /* first vector's dim */
> PetscInt n2=3; /* second vector's dim */
>
> Vec v1; /* n1x1 supervector */
> Vec v2; /* n2x1 reduced-dimension representation of y */
> Mat OP1; /* n1xn2 outer-product matrix */
>
> PetscInt ii; /* for-loop index variables */
>
> PetscInt * locRowIdx1; /* locally-set row indices for VecSetValues */
> PetscScalar * locRowVals1; /* array to hold local vector values */
> PetscInt low1, high1; /* variables to get local-indices' range */
>
> PetscInt * locRowIdx2; /* locally-set row indices for VecSetValues */
> PetscScalar * locRowVals2; /* array to hold local vector values */
> PetscInt low2, high2; /* variables to get local-indices' range */
>
>
> PetscInitialize(&argc,&argv,(char*)0,helpMsg);
>
>
> VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n1,&v1);
> VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n2,&v2);
> VecZeroEntries(v1);
> VecZeroEntries(v2);
>
> MatCreateMPIDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n1,n2,PETSC_NULL,&OP1);
> MatZeroEntries(OP1);
>
>
> VecGetOwnershipRange(v1,&low1,&high1);
> PetscMalloc((high1-low1)*sizeof(PetscInt),&locRowIdx1);
> PetscMalloc((high1-low1)*sizeof(PetscScalar),&locRowVals1);
> /* set locally-set indices and values */
> for (ii = 0; ii < high1-low1; ii++) {
> *(locRowIdx1 + ii) = ii + low1;
> *(locRowVals1 + ii) = pow(2,ii+low1+1);
> } // end of for : set locally-set indices and values
> VecSetValues(v1,high1-low1,locRowIdx1,locRowVals1,INSERT_VALUES);
> PetscFree(locRowIdx1); errCode += PetscFree(locRowVals1);
> /* now, assemble the vector v1 */
> VecAssemblyBegin(v1);
> VecAssemblyEnd(v1);
>
> VecGetOwnershipRange(v2,&low2,&high2);
> PetscMalloc((high2-low2)*sizeof(PetscInt),&locRowIdx2);
> PetscMalloc((high2-low2)*sizeof(PetscScalar),&locRowVals2);
> /* set locally-set indices and values */
> for (ii = 0; ii < high2-low2; ii++) {
> *(locRowIdx2 + ii) = ii + low2;
> *(locRowVals2 + ii) = 1.0 / pow(2,ii+low2);
> } // end of for : set locally-set indices and values
> VecSetValues(v2,high2-low2,locRowIdx2,locRowVals2,INSERT_VALUES);
> PetscFree(locRowIdx2); errCode += PetscFree(locRowVals2);
> /* now, assemble the vector v2 */
> VecAssemblyBegin(v2);
> VecAssemblyEnd(v2);
>
>
> MyMPIVecOuterProd(OP1, v1, v2, INSERT_VALUES);
>
> PetscPrintf(PETSC_COMM_WORLD,"---- vector v1 ----\n");
> VecView(v1,PETSC_VIEWER_STDOUT_WORLD);
> PetscPrintf(PETSC_COMM_WORLD,"---- vector v1 ----\n");
> PetscPrintf(PETSC_COMM_WORLD,"---- vector v2 ----\n");
> VecView(v2,PETSC_VIEWER_STDOUT_WORLD);
> PetscPrintf(PETSC_COMM_WORLD,"---- vector v2 ----\n");
>
> PetscPrintf(PETSC_COMM_WORLD,"---- matrix OP1 ----\n");
> MatView(OP1,PETSC_VIEWER_STDOUT_WORLD);
> PetscPrintf(PETSC_COMM_WORLD,"---- matrix OP1 ----\n");
>
>
> /* destroy OP1 */
> MatDestroy(OP1);
> /* destroy v1 */
> VecDestroy(v1);
> /* destroy v2 */
> VecDestroy(v2);
>
> PetscFinalize();
> return 0;
> }
>
> [hsharma at ifp-32]$ petscmpiexec -np 3 ./OuterProductCheckOutput
>
> ---- vector v1 ----
> Process [0]
> 2
> 4
> Process [1]
> 8
> 16
> Process [2]
> 32
> ---- vector v1 ----
> ---- vector v2 ----
> Process [0]
> 1
> Process [1]
> 0.5
> Process [2]
> 0.25
> ---- vector v2 ----
> ---- matrix OP1 ----
> 2.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
> 4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
> 0.0000000000000000e+00 4.0000000000000000e+00 0.0000000000000000e+00
> 0.0000000000000000e+00 8.0000000000000000e+00 0.0000000000000000e+00
> 0.0000000000000000e+00 0.0000000000000000e+00 8.0000000000000000e+00
> ---- matrix OP1 ----
>
> ---- code and sample output ----
>
More information about the petsc-users
mailing list