[petsc-users] Outer Product of two vectors
Harsh Sharma
hsharma.tgjobs at gmail.com
Wed Oct 6 14:38:44 CDT 2010
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 ----
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101006/c7d918f3/attachment-0001.htm>
More information about the petsc-users
mailing list