Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2779 → Rev 2852

/branches/slepc-3_2-utopia/include/private/matdenseimpl.h
49,6 → 49,8
PetscErrorCode (*setsizes)(MatDense,PetscInt,PetscInt,PetscInt,PetscInt);
PetscErrorCode (*setuppreallocation)(MatDense);
PetscErrorCode (*aresiblings)(MatDense,MatDense,PetscBool*);
PetscErrorCode (*serialize)(MatDense,PetscInt*,PetscScalar*);
PetscErrorCode (*deserialize)(MatDense,PetscInt,PetscScalar*);
PetscErrorCode (*create)(MatDense);
};
 
62,6 → 64,10
PetscBool is_hermitian,is_triangular,is_impl; /* indicates properties of the matrix and how it is stored */
PetscInt matmult_buffersize;
/* suggested size of the operator in bytes for storing the result in MatDenseMatMult_Siblings */
PetscBool use_impl; /* indicates if optimized functions for Hermitian or triangular matrix will be used */
PetscBool use_mpi_queue; /* indicates if the synchronous MPI operations are queued */
PetscInt n_getarray; /* number of unrestored GetArray operations */
PetscBool is_pending; /* indicates if a reduction is been performing */
};
 
/*MC
152,15 → 158,47
} \
}
 
#define PetscCheckMatDenseForRead(A) \
if (!((A)->is_allocated == PETSC_TRUE)) SETERRQ(((PetscObject)(A))->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
PetscErrorCode MatDensePerformReduction(MatDense A,MPI_Op op,PetscBool immediate);
 
#define PetscCheckMatDenseForWrite(A) \
ierr = MatDenseSetUpPreallocation((A));CHKERRQ(ierr); \
if (!((A)->is_allocated == PETSC_TRUE)) SETERRQ(((PetscObject)(A))->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
#undef __FUNCT__
#define __FUNCT__ "PetscCheckMatDenseForRead"
PETSC_STATIC_INLINE PetscErrorCode PetscCheckMatDenseForRead(MatDense A)
{
PetscErrorCode ierr;
PetscFunctionBegin;
if (!(A->is_allocated == PETSC_TRUE)) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
if (A->is_pending && A->n_getarray == 0) {
ierr = MatDensePerformReduction(PETSC_NULL,PETSC_NULL,PETSC_TRUE);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#define PetscCheckMatDenseForUpdate(A) PetscCheckMatDenseForRead(A)
#undef __FUNCT__
#define __FUNCT__ "PetscCheckMatDenseForWrite"
PETSC_STATIC_INLINE PetscErrorCode PetscCheckMatDenseForWrite(MatDense A)
{
PetscErrorCode ierr;
ierr = MatDenseSetUpPreallocation(A);CHKERRQ(ierr);
if (!(A->is_allocated == PETSC_TRUE)) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
if (A->is_pending && A->n_getarray == 0) {
ierr = MatDensePerformReduction(PETSC_NULL,PETSC_NULL,PETSC_TRUE);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PetscCheckMatDenseForUpdate"
PETSC_STATIC_INLINE PetscErrorCode PetscCheckMatDenseForUpdate(MatDense A)
{
PetscErrorCode ierr;
ierr = MatDenseSetUpPreallocation(A);CHKERRQ(ierr);
if (!(A->is_allocated == PETSC_TRUE)) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
if (A->is_pending && A->n_getarray == 0) {
ierr = MatDensePerformReduction(PETSC_NULL,PETSC_NULL,PETSC_TRUE);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#define MatDenseIsHermitian(A) (A->is_hermitian == PETSC_TRUE)
#define MatDenseIsTriangular(A) (A->is_triangular == PETSC_TRUE)
#define MatDenseIsImplicit(A) (A->is_impl == PETSC_TRUE)
/branches/slepc-3_2-utopia/include/slepcmatdense.h
137,7 → 137,10
PetscErrorCode MatDenseGetType(MatDense,const MatDenseType*);
PetscErrorCode MatDenseRegister(const char*,const char*,const char*,PetscErrorCode (*)(MatDense));
PetscErrorCode MatDenseRegisterDestroy(void);
PetscErrorCode MatDenseReduce(MatDense A,MPI_Op op);
 
PetscErrorCode MatDenseCreateBasic(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscScalar *data,MatDense *A);
PetscErrorCode MatDenseCreateBasic(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscScalar *data,MatDense *A);
PetscErrorCode MatDenseSerialize(MatDense A,PetscInt *length,PetscScalar *array);
PetscErrorCode MatDenseDeserialize(MatDense A,PetscInt length,PetscScalar *array);
PETSC_EXTERN_CXX_END
#endif
/branches/slepc-3_2-utopia/src/matdense/interface/basic.c
55,15 → 55,19
PetscFunctionBegin;
PetscValidPointer(A,1);
*A = PETSC_NULL;
ierr = PetscHeaderCreate(B,_p_MatDense,struct _MatDenseOps,MATDENSE_CLASSID,0,"MatDense","Dense Matrix","MatDense",PETSC_COMM_SELF,MatDenseDestroy,MatDenseView);CHKERRQ(ierr);
ierr = PetscHeaderCreate(B,_p_MatDense,struct _MatDenseOps,MATDENSE_CLASSID,0,"MatDense","Dense Matrix","MatDense",comm,MatDenseDestroy,MatDenseView);CHKERRQ(ierr);
*A = B;
B->data = PETSC_NULL;
B->ld = 0;
B->Nmax = B->Nmax = 0;
B->Mmax = B->Nmax = 0;
B->m0 = B->n0 = B->m = B->n = 0;
B->is_allocated = PETSC_FALSE;
B->is_hermitian = B->is_triangular = B->is_impl = PETSC_FALSE;
B->matmult_buffersize = 0;
B->use_impl = PETSC_TRUE;
B->use_mpi_queue = PETSC_FALSE;
B->n_getarray = 0;
B->is_pending = PETSC_FALSE;
PetscFunctionReturn(0);
}
 
119,6 → 123,7
if (n0 < 0 || n0 > A->Nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D invalid or larger than maximum row size %D",n0,A->Nmax);
if (m < 0 || m0+m > A->Mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D invalid or larger than maximum column size %D",m0+m,A->Mmax);
if (n < 0 || n0+n > A->Nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D invalid or larger than maximum row size %D",n0+n,A->Nmax);
if (A->n_getarray > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"The sizes cannot be changed when a MatDenseRestore* call is pending");
if (A->ops->setsizes) {
ierr = (*A->ops->setsizes)(A,m0,n0,m,n);CHKERRQ(ierr);
}
155,8 → 160,8
{
PetscFunctionBegin;
PetscValidHeaderSpecific(mat,MATDENSE_CLASSID,1);
if (m0) *m = mat->m0;
if (n0) *n = mat->n0;
if (m0) *m0 = mat->m0;
if (n0) *n0 = mat->n0;
if (m) *m = mat->m;
if (n) *n = mat->n;
PetscFunctionReturn(0);
185,6 → 190,7
PetscValidHeaderSpecific(*A,MATDENSE_CLASSID,1);
if (--((PetscObject)(*A))->refct > 0) {*A = PETSC_NULL; PetscFunctionReturn(0);}
 
ierr = PetscCheckMatDenseForRead(*A);CHKERRQ(ierr);
if ((*A)->ops->destroy) {
ierr = (*(*A)->ops->destroy)(*A);CHKERRQ(ierr);
}
224,8 → 230,9
PetscValidType(A,1);
PetscValidPointer(v,2);
if (!A->ops->getarray) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
PetscCheckMatDenseForUpdate(A);
ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
ierr = (*A->ops->getarray)(A,v);CHKERRQ(ierr);
A->n_getarray++;
CHKMEMQ;
PetscFunctionReturn(0);
}
256,6 → 263,7
#if defined(PETSC_USE_DEBUG)
CHKMEMQ;
#endif
A->n_getarray--;
if (!A->ops->restorearray) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
ierr = (*A->ops->restorearray)(A,v);CHKERRQ(ierr);
ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr);
296,7 → 304,7
PetscValidType(A,1);
PetscValidPointer(v,2);
if (!A->ops->getarrayread) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
PetscCheckMatDenseForRead(A);
ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);
ierr = (*A->ops->getarrayread)(A,v);CHKERRQ(ierr);
CHKMEMQ;
PetscFunctionReturn(0);
363,7 → 371,7
PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
PetscValidType(A,1);
PetscValidPointer(M,3);
PetscCheckMatDenseForRead(A);
ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);
 
*M = 0;
if (!A->ops->duplicate) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not written for this matrix type");
378,6 → 386,10
B->is_triangular = A->is_triangular;
B->is_impl = A->is_impl;
B->matmult_buffersize = A->matmult_buffersize;
B->use_impl = A->use_impl;
B->use_mpi_queue = A->use_mpi_queue;
B->n_getarray = 0;
B->is_pending = PETSC_FALSE;
}
ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr);
PetscFunctionReturn(0);
455,7 → 467,7
PetscErrorCode ierr;
const char *deft = MATDENSEBASIC;
char type[256];
PetscBool flg;
PetscBool flg,opb;
PetscInt opi;
 
PetscFunctionBegin;
472,6 → 484,10
if (flg) {
A->matmult_buffersize = opi*1024;
}
ierr = PetscOptionsBool("-matdense_enhance_structures","Call structured functions for dense matrices","",A->use_impl,&opb,&flg);CHKERRQ(ierr);
if (flg == PETSC_TRUE) A->use_impl = opb;
ierr = PetscOptionsBool("-matdense_mpi_queue","Queue synchronous MPI calls","",A->use_mpi_queue,&opb,&flg);CHKERRQ(ierr);
if (flg == PETSC_TRUE) A->use_mpi_queue = opb;
 
if (A->ops->setfromoptions) {
ierr = (*A->ops->setfromoptions)(A);CHKERRQ(ierr);
561,6 → 577,7
PetscFunctionReturn(0);
}
 
static PetscErrorCode MatDenseView_Default_ASCII(MatDense xin,PetscViewer viewer);
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseView"
629,7 → 646,7
}
PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
PetscCheckSameComm(mat,1,viewer,2);
if (!mat->is_allocated) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
ierr = PetscCheckMatDenseForRead(mat);CHKERRQ(ierr);
 
ierr = PetscLogEventBegin(MATDENSE_View,mat,viewer,0,0);CHKERRQ(ierr);
ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
646,7 → 663,9
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = (*mat->ops->view)(mat,viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
} else if (!iascii) {
} else if (iascii) {
ierr = MatDenseView_Default_ASCII(mat,viewer);CHKERRQ(ierr);
} else {
SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
}
if (iascii) {
659,7 → 678,67
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseView_Default_ASCII"
static PetscErrorCode MatDenseView_Default_ASCII(MatDense xin,PetscViewer viewer)
{
PetscErrorCode ierr;
PetscInt i,j,n,m,n0,m0;
const char *name;
PetscViewerFormat format;
const PetscScalar *xv;
 
PetscFunctionBegin;
ierr = MatDenseGetSizes(xin,&m0,&n0,&m,&n);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(xin,&xv);CHKERRQ(ierr);
ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
if (format == PETSC_VIEWER_ASCII_MATLAB) {
ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr);
for (i=m0; i<m+m0; i++) {
for (j=n0; j<n+n0; j++) {
PetscScalar v = xv[xin->ld*j+i];
#if defined(PETSC_USE_COMPLEX)
if (PetscImaginaryPart(xv[i]) > 0.0) {
ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",PetscRealPart(v),PetscImaginaryPart(v));CHKERRQ(ierr);
} else if (PetscImaginaryPart(v) < 0.0) {
ierr = PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei ",PetscRealPart(v),-PetscImaginaryPart(v));CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(v));CHKERRQ(ierr);
}
#else
ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double) v);CHKERRQ(ierr);
#endif
}
ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
}
ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr);
} else {
ierr = PetscObjectPrintClassNamePrefixType((PetscObject)xin,viewer,"Vector Object");CHKERRQ(ierr);
for (i=m0; i<m+m0; i++) {
ierr = PetscViewerASCIIPrintf(viewer,"Row %####d: ",i);CHKERRQ(ierr);
for (j=n0; j<n+n0; j++) {
PetscScalar v = xv[xin->ld*j+i];
#if defined(PETSC_USE_COMPLEX)
if (PetscImaginaryPart(v) > 0.0) {
ierr = PetscViewerASCIIPrintf(viewer,"%G + %G i ",PetscRealPart(v),PetscImaginaryPart(v));CHKERRQ(ierr);
} else if (PetscImaginaryPart(v) < 0.0) {
ierr = PetscViewerASCIIPrintf(viewer,"%G - %G i ",PetscRealPart(v),-PetscImaginaryPart(v));CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer,"%G ",PetscRealPart(v));CHKERRQ(ierr);
}
#else
ierr = PetscViewerASCIIPrintf(viewer,"%G ",(double) v);CHKERRQ(ierr);
#endif
}
ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
}
}
ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(xin,&xv);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseFinalizePackage"
/*@C
/branches/slepc-3_2-utopia/src/matdense/interface/ops.c
51,7 → 51,7
PetscErrorCode MatDenseMatMult(MatDense C,PetscScalar alpha,PetscScalar beta,MatDense A,PetscBool At,MatDense B,PetscBool Bt)
{
PetscErrorCode ierr;
PetscInt mA,nA,mB,nB;
PetscInt mA,nA,mB;
PetscBool t,ACsib;
PetscErrorCode (*mult)(MatDense,PetscScalar,PetscScalar,MatDense,PetscBool,MatDense,PetscBool)=PETSC_NULL;
 
65,10 → 65,9
mA = At == PETSC_TRUE ? A->n : A->m;
nA = At == PETSC_TRUE ? A->m : A->n;
mB = Bt == PETSC_TRUE ? B->n : B->m;
nB = Bt == PETSC_TRUE ? B->m : B->n;
PetscCheckMatDenseForRead(A);
PetscCheckMatDenseForRead(B);
PetscCheckMatDenseForWrite(C);
ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);
ierr = PetscCheckMatDenseForRead(B);CHKERRQ(ierr);
ierr = PetscCheckMatDenseForWrite(C);CHKERRQ(ierr);
if (nA!=mB) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",nA,mB);
 
if (mA == 0 || mB == 0) PetscFunctionReturn(0);
100,7 → 99,7
PetscErrorCode MatDenseMatMult_Siblings(MatDense C,PetscScalar alpha,PetscScalar beta,MatDense A,PetscBool At,MatDense B,PetscBool Bt)
{
PetscErrorCode ierr;
PetscInt mA,nA,m0A,n0A,mB,nB,m0C,n0C,mBuff,mBuff0,i;
PetscInt mA,nA,m0A,n0A,nB,m0C,n0C,mBuff,mBuff0,i;
MatDense D;
PetscErrorCode (*mult)(MatDense,PetscScalar,PetscScalar,MatDense,PetscBool,MatDense,PetscBool)=PETSC_NULL;
 
109,7 → 108,6
nA = At == PETSC_TRUE ? A->m : A->n;
m0A = At == PETSC_TRUE ? A->n0 : A->m0;
n0A = At == PETSC_TRUE ? A->m0 : A->n0;
mB = Bt == PETSC_TRUE ? B->n : B->m;
nB = Bt == PETSC_TRUE ? B->m : B->n;
m0C = C->m0;
n0C = C->n0;
151,7 → 149,7
PetscErrorCode MatDenseMatMult_SiblingsA(MatDense C,PetscScalar alpha,PetscScalar beta,MatDense A,PetscBool At,MatDense B,PetscBool Bt)
{
PetscErrorCode ierr;
PetscInt mA,nA,m0A,n0A,mB,nB,m0C,n0C,mBuff,mBuff0,i;
PetscInt mA,nA,m0A,n0A,nB,m0C,n0C,mBuff,mBuff0,i;
MatDense D;
PetscErrorCode (*mult)(MatDense,PetscScalar,PetscScalar,MatDense,PetscBool,MatDense,PetscBool)=PETSC_NULL;
 
160,7 → 158,6
nA = At == PETSC_TRUE ? A->m : A->n;
m0A = At == PETSC_TRUE ? A->n0 : A->m0;
n0A = At == PETSC_TRUE ? A->m0 : A->n0;
mB = Bt == PETSC_TRUE ? B->n : B->m;
nB = Bt == PETSC_TRUE ? B->m : B->n;
m0C = C->m0;
n0C = C->n0;
228,8 → 225,8
PetscValidType(A,1);
PetscValidType(B,2);
PetscCheckSameComm(A,1,B,2);
PetscCheckMatDenseForRead(A);
PetscCheckMatDenseForWrite(B);
ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);
ierr = PetscCheckMatDenseForWrite(B);CHKERRQ(ierr);
if (A->m != B->m || A->n != B->n) SETERRQ4(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: dim (%D,%D) (%D,%D)",A->m,A->n,B->m,B->n);
 
290,8 → 287,8
PetscValidHeaderSpecific(X,MATDENSE_CLASSID,3);
PetscValidHeaderSpecific(Y,MATDENSE_CLASSID,1);
PetscValidLogicalCollectiveScalar(Y,a,2);
PetscCheckMatDenseForRead(X);
PetscCheckMatDenseForUpdate(Y);
ierr = PetscCheckMatDenseForRead(X);CHKERRQ(ierr);
ierr = PetscCheckMatDenseForUpdate(Y);CHKERRQ(ierr);
 
if (X->m != Y->m || X->n != Y->n) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming dense matrix add: %D %D %D %D",X->m,X->n,Y->m,Y->n);
 
426,7 → 423,7
 
PetscFunctionBegin;
PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
PetscCheckMatDenseForUpdate(A);
ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
 
PetscObjectQueryPolymorphicFunction1(A,setexplicit,"MatDenseSetExplicit",&setexplicit,PETSC_TRUE);
ierr = PetscLogEventBegin(MATDENSE_Convert,A,0,0,0);CHKERRQ(ierr);
435,3 → 432,76
PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseSerialize"
/*@
MatDenseSerialize - Copy the matrix elements in a array.
 
Logically Collective on MatDense
 
Input Parameters:
. A - the dense matrix
 
Output Paramenters:
+ lenght - number of scalar to copy or copied (or null)
- array - scalar array with unless length elements (or null)
 
Level: intermediate
 
.keywords: matrix
@*/
PetscErrorCode MatDenseSerialize(MatDense A,PetscInt *length,PetscScalar *array)
{
PetscErrorCode ierr;
PetscErrorCode (*serialize)(MatDense,PetscInt*,PetscScalar*);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
if (length) {
PetscValidPointer(length,2);
}
if (array) {
PetscValidPointer(array,3);
}
ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
 
PetscObjectQueryPolymorphicFunction1(A,serialize,"MatDenseSerialize",&serialize,PETSC_TRUE);
ierr = PetscLogEventBegin(MATDENSE_Convert,A,0,0,0);CHKERRQ(ierr);
ierr = (*serialize)(A,length,array);CHKERRQ(ierr);
ierr = PetscLogEventEnd(MATDENSE_Convert,A,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseDeserialize"
/*@
MatDenseDeserialize - Copy back the elements in the array to the matrix.
 
Logically Collective on MatDense
 
Input Parameters:
+ A - the dense matrix
. lenght - array length (or -1 to avoid checking the size)
- array - scalar array
 
Level: intermediate
 
.keywords: matrix
@*/
PetscErrorCode MatDenseDeserialize(MatDense A,PetscInt length,PetscScalar *array)
{
PetscErrorCode ierr;
PetscErrorCode (*deserialize)(MatDense,PetscInt,PetscScalar*);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
PetscValidPointer(array,3);
ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
 
PetscObjectQueryPolymorphicFunction1(A,deserialize,"MatDenseDeserialize",&deserialize,PETSC_TRUE);
ierr = PetscLogEventBegin(MATDENSE_Convert,A,0,0,0);CHKERRQ(ierr);
ierr = (*deserialize)(A,length,array);CHKERRQ(ierr);
ierr = PetscLogEventEnd(MATDENSE_Convert,A,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/branches/slepc-3_2-utopia/src/matdense/interface/mpi.c New file
0,0 → 1,219
/*
Interface for MatDense class: MPI reduction operations
 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
 
This file is part of SLEPc.
SLEPc is free software: you can redistribute it and/or modify it under the
terms of version 3 of the GNU Lesser General Public License as published by
the Free Software Foundation.
 
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
 
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
 
#include <private/matdenseimpl.h> /*I "slepcmatdense.h" I*/
#include <petscblaslapack.h>
 
/* Structures for MatDenseReduce */
typedef struct MatDenseReductionChunk_ {
MatDense A; /* object matrix to reduce */
PetscScalar *v; /* pointer returned by MatDenseGetArray */
struct MatDenseReductionChunk_ *next; /* next node in this list */
} MatDenseReductionChunk;
 
typedef struct {
MatDenseReductionChunk *r; /* a list of reduction to perform */
PetscInt nr, /* total number of PetscScalar to reduce */
nobjects; /* number of pending objects of this reduction */
MPI_Op op; /* reduction operation */
MPI_Comm comm; /* communicator */
} MatDenseReduction;
 
static PetscErrorCode MatDensePerformReduction_Private(MatDenseReductionChunk *rc0,PetscInt nr,MPI_Op op,MPI_Comm comm,PetscBool free_nr);
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseReduce"
/*@
MatDenseReduce - Performs the MPI_Allreduce with the MPI operation op over a matrix elements.
 
Neighbor-wise Collective on MatDense
 
Input Parameters:
+ A - dense matrix
- op - mpi reduction operation
 
Level: intermediate
 
@*/
PetscErrorCode MatDenseReduce(MatDense A,MPI_Op op)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
if (A->use_mpi_queue) {
ierr = MatDensePerformReduction(A,op,PETSC_FALSE);CHKERRQ(ierr);
} else {
MatDenseReductionChunk rc;
PetscInt l;
ierr = MatDenseSerialize(A,&l,PETSC_NULL);CHKERRQ(ierr);
ierr = MatDenseGetArray(A,&rc.v);CHKERRQ(ierr);
rc.A = A;
rc.next = PETSC_NULL;
ierr = MatDensePerformReduction_Private(&rc,l,op,((PetscObject)A)->comm,PETSC_FALSE);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#include <omp.h>
 
#undef __FUNCT__
#define __FUNCT__ "MatDensePerformReduction"
PetscErrorCode MatDensePerformReduction(MatDense A,MPI_Op op,PetscBool immediate)
{
PetscErrorCode ierr;
PetscInt l;
MatDenseReductionChunk nr0,*nr;
PetscScalar *v;
static MatDenseReduction *currentr = PETSC_NULL; /* reduction queue */
static PetscBool locks_initialized = PETSC_FALSE;
static omp_lock_t lock; /* reduction in curse */
static omp_lock_t lockq; /* processing the queue */
PetscBool new_thread,waitq;
 
PetscFunctionBegin;
/* Quick return */
if (!A && (!locks_initialized || !immediate)) PetscFunctionReturn(0);
 
/* Init lock if needed */
if (!locks_initialized) {
omp_init_lock(&lock);
omp_init_lock(&lockq);
locks_initialized = PETSC_TRUE;
}
 
if (A) {
/* Prepare reduction */
ierr = MatDenseSerialize(A,&l,PETSC_NULL);CHKERRQ(ierr);
ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
A->is_pending = PETSC_TRUE;
nr0.A = A;
nr0.v = v;
nr0.next = PETSC_NULL;
if (omp_test_lock(&lock)) {
ierr = MatDensePerformReduction_Private(&nr0,l,op,((PetscObject)A)->comm,PETSC_FALSE);CHKERRQ(ierr);
omp_unset_lock(&lock);
PetscFunctionReturn(0);
}
#pragma omp critical lock_currentr
{
/* If the communicator or the operation of the currect reduction are not compatible the reduction is not accumulated to them */
waitq = (currentr && (((PetscObject)A)->comm != currentr->comm || currentr->op != op))?PETSC_TRUE:PETSC_FALSE;
}
if (waitq) {
omp_set_lock(&lockq);
omp_unset_lock(&lockq);
}
#pragma omp critical lock_currentr
{
/* Create the list of reductions if needed */
if (!currentr) {
ierr = PetscMalloc(sizeof(MatDenseReduction),&currentr);CHKERRQ(ierr);
currentr->r = PETSC_NULL;
currentr->nr = 0;
currentr->nobjects = 0;
currentr->op = op;
currentr->comm = ((PetscObject)A)->comm;
new_thread = PETSC_TRUE;
} else {
new_thread = PETSC_FALSE;
}
/* Accumulate reduction */
ierr = PetscMalloc(sizeof(MatDenseReductionChunk),&nr);CHKERRQ(ierr);
nr->A = A;
nr->v = v;
nr->next = currentr->r;
currentr->r = nr;
currentr->nobjects++;
currentr->nr+= l;
}
/* Perform the reduction if it is immediate */
if (new_thread) {
omp_set_lock(&lockq);
#pragma omp parallel
#pragma omp single nowait
{
PetscErrorCode ierr;
MatDenseReduction *currentr0;
omp_set_lock(&lock);
#pragma omp critical lock_currentr
{
currentr0 = currentr;
currentr = PETSC_NULL;
}
omp_unset_lock(&lockq);
ierr = MatDensePerformReduction_Private(currentr0->r,currentr0->nr,currentr0->op,currentr0->comm,PETSC_TRUE);CHKERRABORT(currentr0->comm,ierr);
ierr = PetscFree(currentr0);CHKERRABORT(currentr0->comm,ierr);
omp_unset_lock(&lock);
}
}
}
 
if (immediate) {
omp_set_lock(&lockq);
omp_unset_lock(&lockq);
omp_set_lock(&lock);
omp_unset_lock(&lock);
}
 
PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "MatDensePerformReduction_Private"
static PetscErrorCode MatDensePerformReduction_Private(MatDenseReductionChunk *rc0,PetscInt nr,MPI_Op op,MPI_Comm comm,PetscBool free_nr)
{
PetscErrorCode ierr;
PetscInt l;
MatDenseReductionChunk *oldr,*rc;
PetscScalar *in,*v;
 
PetscFunctionBegin;
 
ierr = PetscMalloc(sizeof(PetscScalar)*nr,&in);CHKERRQ(ierr);
/* Copy the matrix elements into the array */
for(rc = rc0, v = in; rc; rc = rc->next, v+= l) {
ierr = MatDenseSerialize(rc->A,&l,v);CHKERRQ(ierr);
}
/* Perform the reduction */
MPI_Allreduce(MPI_IN_PLACE,in,nr,MPIU_SCALAR,op,comm);
/* Copy back the reduced elements and destroy the list */
for(rc = rc0, v = in; rc; ) {
ierr = MatDenseDeserialize(rc->A,-1,v);CHKERRQ(ierr);
rc->A->is_pending = PETSC_FALSE;
ierr = MatDenseRestoreArray(rc->A,&rc->v);CHKERRQ(ierr);
oldr = rc;
rc = rc->next;
if (free_nr) {
ierr = PetscFree(oldr);CHKERRQ(ierr);
}
}
ierr = PetscFree(in);CHKERRQ(ierr);
 
PetscFunctionReturn(0);
}
/branches/slepc-3_2-utopia/src/matdense/interface/makefile
23,10 → 23,10
 
CFLAGS =
FFLAGS =
SOURCEC = basic.c itregis.c ops.c
SOURCEC = basic.c itregis.c ops.c mpi.c
SOURCEF =
SOURCEH =
OBJSC = basic.o itregis.o ops.o
OBJSC = basic.o itregis.o ops.o mpi.o
LIBBASE = libslepc
DIRS =
MANSEC = MatDense
/branches/slepc-3_2-utopia/src/matdense/impls/basic/basic.c
26,8 → 26,7
 
typedef struct {
PetscScalar *v;
PetscBool use_impl,
user_allocated;
PetscBool user_allocated;
PetscInt siblings;
} MatDense_Basic;
 
35,19 → 34,134
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseBlasMatMult_Basic"
/*
A ------C B ------C ret C --- Operations
ExplHermN ExplHermN ExplNormN hemm L | hemm R
ExplHermN ExplHermC "" hemm R
ExplHermN ImplHermN "" hemm R
ExplHermN ImplHermC "" hemm R
ExplHermN ExplTriaN "" hemm L | trmm R
ExplHermN ExplTriaC "" trmm R
ExplHermN ImplTriaN "" trmm R
ExplHermN ImplTriaC "" trmm R
ExplHermN ExplNormN "" hemm L
ExplHermN ExplNormC "" gemm
 
ExplHermC ExplHermN "" hemm L
ExplHermC ExplHermC "" gemm
ExplHermC ImplHermN "" expl B, gemm
ExplHermC ImplHermC "" expl B, gemm
ExplHermC ExplTriaN "" hemm L
ExplHermC ExplTriaC "" gemm
ExplHermC ImplTriaN "" expl B, gemm
ExplHermC ImplTriaC "" expl B, gemm
ExplHermC ExplNormN "" hemm L
ExplHermC ExplNormC "" gemm
 
ImplHermN ExplHermN "" hemm L
ImplHermN ExplHermC "" expl A, hemm R | expl A, gemm
ImplHermN ImplHermN "" expl B, hemm L | expl AB, gemm
ImplHermN ImplHermC "" expl A, hemm R | expl AB, gemm
ImplHermN ExplTriaN "" hemm L
ImplHermN ExplTriaC "" expl A, trmm R | expl A, gemm
ImplHermN ImplTriaN "" expl B, hemm L | expl A, trmm R | expl AB, gemm
ImplHermN ImplTriaC "" expl A, trmm R | expl AB, gemm
ImplHermN ExplNormN "" hemm L
ImplHermN ExplNormC "" expl A, gemm
 
ImplHermC ExplHermN "" hemm L
ImplHermC ExplHermC "" expl A, hemm R | expl A, gemm
ImplHermC ImplHermN "" expl B, hemm L | expl AB, gemm
ImplHermC ImplHermC "" expl AB, gemm
ImplHermC ExplTriaN "" hemm L
ImplHermC ExplTriaC "" expl A, gemm
ImplHermC ImplTriaN "" expl B, hemm L | expl AB, gemm
ImplHermC ImplTriaC "" expl AB, gemm
ImplHermC ExplNormN "" hemm L
ImplHermC ExplNormC "" expl A, gemm
 
ExplTriaN ExplHermN "" hemm R | trmm L
ExplTriaN ExplHermC "" hemm R
ExplTriaN ImplHermN "" hemm R
ExplTriaN ImplHermC "" hemm R
ExplTriaN ExplTriaN "" trmm L | trmm R
ExplTriaN ExplTriaC "" trmm R
ExplTriaN ImplTriaN "" trmm R
ExplTriaN ImplTriaC "" trmm R
ExplTriaN ExplNormN "" trmm L
ExplTriaN ExplNormC "" trmm R
 
ExplTriaC ExplHermN "" trmm L
ExplTriaC ExplHermC "" gemm
ExplTriaC ImplHermN "" expl B, trmm L | expl B, gemm
ExplTriaC ImplHermC "" expl B, gemm
ExplTriaC ExplTriaN "" trmm L
ExplTriaC ExplTriaC "" gemm
ExplTriaC ImplTriaN "" trmm L
ExplTriaC ImplTriaC "" expl B, gemm
ExplTriaC ExplNormN "" trmm L
ExplTriaC ExplNormC "" gemm
 
ImplTriaN ExplHermN "" trmm L
ImplTriaN ExplHermC "" expl A, hemm R | expl A, gemm
ImplTriaN ImplHermN "" expl B, trmm L | expl A, hemm R | expl AB, gemm
ImplTriaN ImplHermC "" expl A, hemm R | expl AB, gemm
ImplTriaN ExplTriaN "" trmm L
ImplTriaN ExplTriaC "" expl A, trmm R | expl A, gemm
ImplTriaN ImplTriaN "" expl B, trmm L | expl AB, gemm
ImplTriaN ImplTriaC "" expl A, trmm R | expl AB, gemm
ImplTriaN ExplNormN "" trmm L
ImplTriaN ExplNormC "" expl A, gemm
 
ImplTriaC ExplHermN "" trmm L
ImplTriaC ExplHermC "" expl A, gemm
ImplTriaC ImplHermN "" expl A, hemm R | expl AB, gemm
ImplTriaC ImplHermC "" expl AB, gemm
ImplTriaC ExplTriaN "" trmm L
ImplTriaC ExplTriaC "" expl A, gemm
ImplTriaC ImplTriaN "" expl B, trmm L | expl AB, gemm
ImplTriaC ImplTriaC "" expl AB, gemm
ImplTriaC ExplNormN "" trmm L
ImplTriaC ExplNormC "" expl A, gemm
 
ExplNormN ExplHermN "" hemm R
ExplNormN ExplHermC "" hemm R
ExplNormN ImplHermN "" hemm R
ExplNormN ImplHermC "" hemm R
ExplNormN ExplTriaN "" trmm R
ExplNormN ExplTriaC "" trmm R
ExplNormN ImplTriaN "" trmm R
ExplNormN ImplTriaC "" trmm R
ExplNormN ExplNormN "" gemm
ExplNormN ExplNormC "" gemm
 
ExplNormC ExplHermN "" gemm
ExplNormC ExplHermC "" gemm
ExplNormC ImplHermN "" expl B, gemm
ExplNormC ImplHermC "" expl B, gemm
ExplNormC ExplTriaN "" gemm
ExplNormC ExplTriaC "" gemm
ExplNormC ImplTriaN "" expl B, gemm
ExplNormC ImplTriaC "" expl B, gemm
ExplNormC ExplNormN "" gemm
ExplNormC ExplNormC "" gemm
*/
PetscErrorCode MatDenseBlasMatMult_Basic(MatDense C,PetscScalar alpha,PetscScalar beta,MatDense A,PetscBool At,MatDense B,PetscBool Bt)
{
PetscErrorCode ierr;
PetscBLASInt mA,nA,mB,nB,ldA,ldB,ldC;
PetscScalar *vC;
PetscScalar *vC_,*vC;
const PetscScalar
*vA, *vB;
*vA,*vA_,*vB,*vB_;
const char *N = "N", *T = "C", *qA = N, *qB = N;
 
PetscFunctionBegin;
ierr = MatDenseGetArrayRead(A,&vA);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(B,&vB);CHKERRQ(ierr);
ierr = MatDenseGetArray(C,&vC);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(A,&vA_);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(B,&vB_);CHKERRQ(ierr);
ierr = MatDenseGetArray(C,&vC_);CHKERRQ(ierr);
vA = &vA_[A->ld*A->n0+A->m0];
vB = &vB_[B->ld*B->n0+B->m0];
vC = &vC_[C->ld*C->n0+C->m0];
 
mA = At == PETSC_TRUE ? A->n : A->m;
nA = At == PETSC_TRUE ? A->m : A->n;
65,25 → 179,27
else if (!At && Bt) *vC = *vA * PetscConj(*vB);
else if (At && Bt) *vC = PetscConj(*vA) * PetscConj(*vB);
ierr = PetscLogFlops(1);CHKERRQ(ierr);
} else if (!A->is_impl && !B->is_impl) {
BLASgemm_(qA,qB,&mA,&nB,&nA,&alpha,(PetscScalar*)vA,&ldA,(PetscScalar*)vB,&ldB,&beta,vC,&ldC);
ierr = PetscLogFlops(mA*nB*2.0*nA);CHKERRQ(ierr);
/* A Hermitian, B !impl !trans) */
} else if (A->is_hermitian && !B->is_impl && !Bt) {
/* A Hermitian, B !trans) */
} else if (A->use_impl && A->is_hermitian && !Bt) {
if (B->is_impl) { ierr = MatDenseSetExplicit(B);CHKERRQ(ierr); }
BLAShemm_("L","U",&mA,&nB,&alpha,(PetscScalar*)vA,&ldA,(PetscScalar*)vB,&ldB,&beta,vC,&ldC);
ierr = PetscLogFlops(mA*nB*nA);CHKERRQ(ierr);
/* A !impl !trans, B Hermitian) */
} else if (!A->is_impl && !At && B->is_hermitian) {
/* A !trans, B Hermitian) */
} else if (A->use_impl && !At && B->is_hermitian) {
if (A->is_impl) { ierr = MatDenseSetExplicit(A);CHKERRQ(ierr); }
BLAShemm_("R","U",&mA,&nB,&alpha,(PetscScalar*)vB,&ldB,(PetscScalar*)vA,&ldA,&beta,vC,&ldC);
ierr = PetscLogFlops(mA*nB*nA);CHKERRQ(ierr);
/* A triangular, B !impl !trans) */
} else if (A->is_triangular && !B->is_impl && !Bt) {
/* A triangular, B !trans) */
} else if (A->use_impl && A->is_triangular && !Bt) {
if (B->is_impl) { ierr = MatDenseSetExplicit(B);CHKERRQ(ierr); }
BLAStrmm_("L","U",qA,"N",&mB,&nB,&alpha,(PetscScalar*)vA,&ldA,(PetscScalar*)vB,&ldB);
ierr = PetscLogFlops(mA*nB*nA);CHKERRQ(ierr);
/* A !impl !trans, B triangular) */
} else if (!A->is_impl && !At && B->is_triangular) {
/* A !trans, B triangular) */
} else if (A->use_impl && !At && B->is_triangular) {
if (A->is_impl) { ierr = MatDenseSetExplicit(A);CHKERRQ(ierr); }
BLAStrmm_("R","U",qB,"N",&mA,&nA,&alpha,(PetscScalar*)vB,&ldB,(PetscScalar*)vA,&ldA);
ierr = PetscLogFlops(mA*nB*nA);CHKERRQ(ierr);
/* otherwise */
} else {
if (A->is_impl) { ierr = MatDenseSetExplicit(A);CHKERRQ(ierr); }
if (B->is_impl) { ierr = MatDenseSetExplicit(B);CHKERRQ(ierr); }
92,9 → 208,9
}
C->is_hermitian = C->is_triangular = C->is_impl = PETSC_FALSE;
 
ierr = MatDenseRestoreArray(C,&vC);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(B,&vB);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(A,&vA);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(C,&vC_);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(B,&vB_);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(A,&vA_);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
103,17 → 219,19
PetscErrorCode MatDenseBlasCopy_Basic(PetscInt n,MatDense A,PetscInt dA,PetscInt iA,MatDense B,PetscInt dB,PetscInt iB)
{
PetscErrorCode ierr;
PetscScalar *vB;
PetscScalar *vB,*vB_;
const PetscScalar
*vA;
*vA,*vA_;
PetscBLASInt n_=PetscBLASIntCast(n),incA=PetscBLASIntCast(iA),incB=PetscBLASIntCast(iB);
 
PetscFunctionBegin;
ierr = MatDenseGetArrayRead(A,&vA);CHKERRQ(ierr);
ierr = MatDenseGetArray(B,&vB);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(A,&vA_);CHKERRQ(ierr);
ierr = MatDenseGetArray(B,&vB_);CHKERRQ(ierr);
vA = &vA_[A->ld*A->n0+A->m0];
vB = &vB_[B->ld*B->n0+B->m0];
BLAScopy_(&n_,(PetscScalar*)&vA[dA],&incA,&vB[dB],&incB);
ierr = MatDenseRestoreArray(B,&vB);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(A,&vA);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(B,&vB_);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(A,&vA_);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
122,17 → 240,19
PetscErrorCode MatDenseBlasAXPY_Basic(PetscInt n,PetscScalar a,MatDense X,PetscInt dX,PetscInt iX,MatDense Y,PetscInt dY,PetscInt iY)
{
PetscErrorCode ierr;
PetscScalar *vY;
PetscScalar *vY,*vY_;
const PetscScalar
*vX;
*vX,*vX_;
PetscBLASInt n_=PetscBLASIntCast(n),incX=PetscBLASIntCast(iX),incY=PetscBLASIntCast(iY);
 
PetscFunctionBegin;
ierr = MatDenseGetArrayRead(X,&vX);CHKERRQ(ierr);
ierr = MatDenseGetArray(Y,&vY);CHKERRQ(ierr);
ierr = MatDenseGetArrayRead(X,&vX_);CHKERRQ(ierr);
ierr = MatDenseGetArray(Y,&vY_);CHKERRQ(ierr);
vX = &vX_[X->ld*X->n0+X->m0];
vY = &vY_[Y->ld*Y->n0+Y->m0];
BLASaxpy_(&n_,&a,(PetscScalar*)&vX[dX],&incX,&vY[dY],&incY);
ierr = MatDenseRestoreArray(Y,&vY);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(X,&vX);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(Y,&vY_);CHKERRQ(ierr);
ierr = MatDenseRestoreArrayRead(X,&vX_);CHKERRQ(ierr);
ierr = PetscLogFlops(2*n);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
142,12 → 262,13
PetscErrorCode MatDenseSetExplicit_Basic(MatDense A)
{
PetscErrorCode ierr;
PetscScalar *vA;
PetscScalar *vA,*vA_;
PetscInt i;
PetscBLASInt n,incX,one=1;
 
PetscFunctionBegin;
ierr = MatDenseGetArray(A,&vA);CHKERRQ(ierr);
ierr = MatDenseGetArray(A,&vA_);CHKERRQ(ierr);
vA = &vA_[A->ld*A->n0+A->m0];
if (MatDenseIsImplicitHermitian(A)) {
for (i=0; i<A->n; i++) {
n=A->n-i-1; incX=A->ld;
160,7 → 281,7
}
A->is_impl = PETSC_FALSE;
 
ierr = MatDenseRestoreArray(A,&vA);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(A,&vA_);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
187,7 → 308,6
#define __FUNCT__ "MatDenseGetArray_Basic"
PetscErrorCode MatDenseGetArray_Basic(MatDense A,PetscScalar *v[])
{
PetscErrorCode ierr;
MatDense_Basic *data = (MatDense_Basic*)A->data;
 
PetscFunctionBegin;
199,8 → 319,6
#define __FUNCT__ "MatDenseRestoreArray_Basic"
PetscErrorCode MatDenseRestoreArray_Basic(MatDense A,PetscScalar *v[])
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscFunctionReturn(0);
}
209,7 → 327,6
#define __FUNCT__ "MatDenseGetArrayRead_Basic"
PetscErrorCode MatDenseGetArrayRead_Basic(MatDense A,const PetscScalar *v[])
{
PetscErrorCode ierr;
MatDense_Basic *data = (MatDense_Basic*)A->data;
 
PetscFunctionBegin;
221,8 → 338,6
#define __FUNCT__ "MatDenseRestoreArrayRead_Basic"
PetscErrorCode MatDenseRestoreArrayRead_Basic(MatDense A,const PetscScalar *v[])
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscFunctionReturn(0);
}
240,12 → 355,10
ierr = MatDenseSetMaxSizes(*M,A->m,A->n);CHKERRQ(ierr);
ierr = MatDenseSetSizes(*M,0,0,A->m,A->n);CHKERRQ(ierr);
ierr = MatDenseCreateBasic_Private(*M,op == MATDENSE_MAKE_SIBLING ? data : PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
((MatDense_Basic*)(*M)->data)->use_impl = data->use_impl;
if (op == MATDENSE_COPY_VALUES) {
ierr = MatDenseSetUpPreallocation(*M);CHKERRQ(ierr);
ierr = MatDenseCopy(A,*M);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
254,15 → 367,7
#define __FUNCT__ "MatDenseSetFromOptions_Basic"
PetscErrorCode MatDenseSetFromOptions_Basic(MatDense A)
{
PetscErrorCode ierr;
PetscBool flg,set;
MatDense_Basic *data = (MatDense_Basic*)A->data;
 
PetscFunctionBegin;
 
ierr = PetscOptionsBool("-matdense_enhance_structures","Call structured functions for dense matrices","",data->use_impl,&set,&flg);CHKERRQ(ierr);
if (flg == PETSC_TRUE) data->use_impl = set;
 
PetscFunctionReturn(0);
}
 
281,6 → 386,72
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseSerialize_Basic"
PetscErrorCode MatDenseSerialize_Basic(MatDense A,PetscInt *length,PetscScalar *array)
{
PetscErrorCode ierr;
const PetscScalar *vA,*vA_;
PetscInt i;
PetscBLASInt n,one=1;
 
PetscFunctionBegin;
ierr = MatDenseGetArrayRead(A,&vA_);CHKERRQ(ierr);
vA = &vA_[A->ld*A->n0+A->m0];
if (A->use_impl && A->m == A->n && (A->is_hermitian || A->is_triangular)) {
if (length) *length = A->n*(A->n+1)/2;
if (array) {
for (i=0; i<A->n; i++) {
n = i+1;
BLAScopy_(&n,&vA[A->ld*i],&one,array,&one);
array+= i+1;
}
}
} else {
if (length) *length = A->n*A->m;
if (array) {
n = PetscBLASIntCast(A->m);
for (i=0; i<A->n; i++) {
BLAScopy_(&n,&vA[A->ld*i],&one,array,&one);
array+= A->m;
}
}
}
ierr = MatDenseRestoreArrayRead(A,&vA_);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatDenseDeserialize_Basic"
PetscErrorCode MatDenseDeserialize_Basic(MatDense A,PetscInt length,PetscScalar *array)
{
PetscErrorCode ierr;
PetscScalar *vA,*vA_;
PetscInt i;
PetscBLASInt n,one=1;
 
PetscFunctionBegin;
ierr = MatDenseGetArray(A,&vA_);CHKERRQ(ierr);
vA = &vA_[A->ld*A->n0+A->m0];
if (A->use_impl && A->m == A->n && (A->is_hermitian || A->is_triangular)) {
if (length >= 0 && length != A->n*(A->n+1)/2) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Invalid array length to deseralize, that has to be %D!",A->n*(A->n+1)/2);
for (i=0; i<A->n; i++) {
n = i+1;
BLAScopy_(&n,array,&one,&vA[A->ld*i],&one);
array+= i+1;
}
} else {
if (length >= 0 && length != A->n*A->m) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Invalid array length to deseralize, that has to be %D!",A->n*A->m);
n = PetscBLASIntCast(A->m);
for (i=0; i<A->n; i++) {
BLAScopy_(&n,array,&one,&vA[A->ld*i],&one);
array+= A->m;
}
}
ierr = MatDenseRestoreArray(A,&vA_);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
static struct _MatDenseOps MatDenseOps_Basic = {
/* 0*/
MatDenseGetArray_Basic,
304,6 → 475,8
0/*MatDenseSetSizes_Basic*/,
MatDenseSetUpPreallocation_Basic,
0/*MatDenseAreSiblings_Basic*/,
MatDenseSerialize_Basic,
MatDenseDeserialize_Basic,
0
};
 
319,7 → 492,6
ierr = PetscNewLog(A,MatDense_Basic,&b);CHKERRQ(ierr);
b->v = v;
b->user_allocated = v?PETSC_TRUE:PETSC_FALSE;
b->use_impl = PETSC_FALSE;
b->siblings = 1;
}
A->data = (void*)b;
/branches/slepc-3_2-utopia/src/matdense/examples/tests/test2.c New file
0,0 → 1,97
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
 
This file is part of SLEPc.
SLEPc is free software: you can redistribute it and/or modify it under the
terms of version 3 of the GNU Lesser General Public License as published by
the Free Software Foundation.
 
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
 
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
 
static char help[] = "Tests MatDenseReduce.\n\n";
#include <slepcmatdense.h>
 
#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv )
{
MatDense A,B; /* matrices */
PetscScalar *v;
PetscInt n=45,m,i,j;
PetscBool flag,verbose = PETSC_FALSE;
PetscErrorCode ierr;
 
SlepcInitialize(&argc,&argv,(char*)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);CHKERRQ(ierr);
if(!flag) m=n;
ierr = PetscOptionsGetBool(PETSC_NULL,"-verbose",&verbose,&flag);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nDense matrices of size %Dx%D\n\n",n,m);CHKERRQ(ierr);
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create matrices
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatDenseCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatDenseSetMaxSizes(A,m,n);CHKERRQ(ierr);
ierr = MatDenseSetSizes(A,0,0,m,n);CHKERRQ(ierr);
ierr = MatDenseSetFromOptions(A);CHKERRQ(ierr);
ierr = MatDenseSetUpPreallocation(A);CHKERRQ(ierr);
ierr = MatDenseDuplicate(A,MATDENSE_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr);
ierr = MatDenseSetFromOptions(B);CHKERRQ(ierr);
ierr = MatDenseSetUpPreallocation(B);CHKERRQ(ierr);
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Initialize matrices
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
for (j=0; j<n; j++) {
for (i=0; i<m; i++) {
v[m*j+i] = 1.0;
}
}
ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
if (verbose) {
ierr = MatDenseView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
 
ierr = MatDenseGetArray(B,&v);CHKERRQ(ierr);
for (j=0; j<n; j++) {
for (i=0; i<m; i++) {
v[m*j+i] = i % 2 ? 1.0 : 2.0;
}
}
ierr = MatDenseRestoreArray(B,&v);CHKERRQ(ierr);
if (verbose) {
ierr = MatDenseView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Test operations
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatDenseReduce(A,MPIU_SUM);CHKERRQ(ierr);
ierr = MatDenseReduce(B,MPIU_SUM);CHKERRQ(ierr);
ierr = MatDenseAXPY(B,1.0,A);CHKERRQ(ierr);
if (verbose) {
ierr = MatDenseView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
 
ierr = MatDenseDestroy(&A);CHKERRQ(ierr);
ierr = MatDenseDestroy(&B);CHKERRQ(ierr);
ierr = SlepcFinalize();CHKERRQ(ierr);
return 0;
}
 
/branches/slepc-3_2-utopia/src/matdense/examples/tests/makefile
27,9 → 27,10
EXAMPLESC = test1.c
EXAMPLESF =
MANSEC = MatDense
TESTS = test1
TESTS = test1 test2
 
TESTEXAMPLES_C = test1.PETSc runtest1_1 test1.rm
TESTEXAMPLES_C = test1.PETSc runtest1_1 test1.rm \
test2.PETSc runtest2_1 test2.rm
 
TESTEXAMPLES_BLOPEX =
 
39,6 → 40,10
-${CLINKER} -o test1 test1.o ${SLEPC_LIB}
${RM} test1.o
 
test2: test2.o chkopts
-${CLINKER} -o test2 test2.o ${SLEPC_LIB} -fopenmp
${RM} test2.o
 
#------------------------------------------------------------------------------------
runtest1_1:
-@touch test1_1.tmp; \
47,3 → 52,16
else echo "Possible problem with test1_1, diffs above"; fi; \
${RM} -f test1_1.tmp
 
runtest2_1:
-@touch test2_1.tmp; \
${MPIEXEC} -np 1 ./test2 > test2_1.tmp 2>&1; \
if (${DIFF} output/test2_1.out test2_1.tmp) then true; \
else echo "Possible problem with test2_1, diffs above"; fi; \
${RM} -f test2_1.tmp
 
runtest2_2:
-@touch test2_2.tmp; \
${MPIEXEC} -np 2 ./test2 > test2_2.tmp 2>&1; \
if (${DIFF} output/test2_2.out test2_2.tmp) then true; \
else echo "Possible problem with test2_2, diffs above"; fi; \
${RM} -f test2_2.tmp