Subversion Repositories slepc-dev

Rev

Blame | Compare with Previous | Last modification | View Log | RSS feed

/*
   Interface for MatDense class.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   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>

PetscFList       MatDenseList = 0;
PetscBool        MatDenseRegisterAllCalled = PETSC_FALSE;
PetscClassId     MATDENSE_CLASSID = 0;
PetscLogEvent    MATDENSE_Duplicate = 0,MATDENSE_SetUpPreallocation = 0,MATDENSE_MatMult = 0,MATDENSE_Copy = 0,MATDENSE_AXPY = 0,MATDENSE_View = 0,MATDENSE_Convert = 0,MATDENSE_BlasMatMult = 0;
static PetscBool MatDensePackageInitialized = PETSC_FALSE;

#undef __FUNCT__
#define __FUNCT__ "MatDenseCreate"
/*@
   MatDenseCreate - Creates an empty dense matrix object. The type can then be
   set with MatDenseSetType(), or MatDenseSetFromOptions().

   Collective on MatDense

   Input Parameters:
.  m - the dense matrix

   Note:
   If you never call MatDenseSetType() or MatDenseFromOptions() it will generate
   an error when you try to use the object.

   Level: developer
@*/

PetscErrorCode MatDenseCreate(MPI_Comm comm,MatDense *A)
{
  MatDense        B;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  PetscValidPointer(A,1);
  *A = PETSC_NULL;
  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->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);
}

#undef __FUNCT__
#define __FUNCT__ "MatDenseSetMaxSizes"
/*@
   MatDenseSetMaxSizes - Sets the maximum size of rows (also called
   leading dimension) and columns.

   Collective on MatDense

   Input Parameters:
+  A - the dense matrix
.  M - the maximum rows size
-  N - the maximum columns size

   Level: developer
@*/

PetscErrorCode MatDenseSetMaxSizes(MatDense A, PetscInt M, PetscInt N)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  if (A->is_allocated == PETSC_TRUE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset sizes of the matrix after set the data");
  if (M <= 0 || N <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid dimension value");
  A->ld = A->Mmax = M;
  A->Nmax = N;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseSetSizes"
/*@
  MatDenseSetSizes - Sets the local sizes.

  Collective on MatDense

  Input Parameters:
+  A - the matrix
.  m0 - displacement of the first row
.  n0 - displacement of the first column
.  m - number of rows
-  n - number of columns

  Level: beginner
@*/

PetscErrorCode MatDenseSetSizes(MatDense A,PetscInt m0,PetscInt n0,PetscInt m,PetscInt n)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  if (m0 < 0 || m0 > A->Mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D invalid or larger than maximum column size %D",m0,A->Mmax);
  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);
  }
  A->m = m;
  A->n = n;
  A->m0 = m0;
  A->n0 = n0;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseGetSizes"
/*@
   MatDenseGetSizes - Returns the current numbers of rows and columns in a matrix.

   Not Collective

   Input Parameter:
.  mat - the matrix

   Output Parameters:
+  m0 - displacement of the first row
.  n0 - displacement of the first column
.  m - the number of rows
-  n - the number of columns

   Note: both output parameters can be PETSC_NULL on input.

   Level: beginner

   Concepts: matrices^size
@*/

PetscErrorCode  MatDenseGetSizes(MatDense mat,PetscInt *m0,PetscInt *n0,PetscInt *m,PetscInt* n)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MATDENSE_CLASSID,1);
  if (m0) *m0 = mat->m0;
  if (n0) *n0 = mat->n0;
  if (m) *m = mat->m;
  if (n) *n = mat->n;
  PetscFunctionReturn(0);
}


#undef __FUNCT__
#define __FUNCT__ "MatDenseDestroy"
/*@
   MatDenseDestroy - Frees space taken by a matrix.

   Collective on MatDense

   Input Parameter:
.  A - the matrix

   Level: beginner

@*/

PetscErrorCode  MatDenseDestroy(MatDense *A)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (!*A) PetscFunctionReturn(0);
  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);
  }

  ierr = PetscHeaderDestroy(A);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseGetArray"
/*@C
   MatDenseGetArray - Returns a pointer to the element values in the matrix.
   The result of this routine is dependent on the underlying matrix data
   structure, and may not even work for certain matrix types.  You MUST
   call MatRestoreArray() when you no longer need to access the array.

   Not Collective

   Input Parameter:
.  A - the matrix

   Output Parameter:
.  v - the location of the values

   Level: advanced

   Concepts: matrices^access array

.seealso: MatDenseRestoreArray()
@*/

PetscErrorCode  MatDenseGetArray(MatDense A,PetscScalar *v[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidPointer(v,2);
  if (!A->ops->getarray) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
  ierr = PetscCheckMatDenseForUpdate(A);CHKERRQ(ierr);
  ierr = (*A->ops->getarray)(A,v);CHKERRQ(ierr);
  A->n_getarray++;
  CHKMEMQ;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseRestoreArray"
/*@C
   MatDenseRestoreArray - Restores the matrix after MatDenseGetArray() has been called.

   Not Collective

   Input Parameter:
+  A - the matrix
-  v - the location of the values

   Level: advanced

.seealso: MatGetArray(), MatRestoreArrayF90()
@*/

PetscErrorCode  MatDenseRestoreArray(MatDense A,PetscScalar *v[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidPointer(v,2);
#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);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseGetArrayRead"
/*@C
   MatDenseGetArrayRead - Returns a pointer to the element values in the matrix.
   The result of this routine is dependent on the underlying matrix data
   structure, and may not even work for certain matrix types.  You MUST
   call MatRestoreArrayRead() when you no longer need to access the array.

   Not Collective

   Input Parameter:
.  A - the matrix

   Output Parameter:
.  v - the location of the values

   Level: advanced

   Note:
   The routine considers that the values of the array will not be modified.

   Concepts: matrices^access array

.seealso: MatDenseRestoreArrayRead()
@*/

PetscErrorCode  MatDenseGetArrayRead(MatDense A,const PetscScalar *v[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidPointer(v,2);
  if (!A->ops->getarrayread) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
  ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);
  ierr = (*A->ops->getarrayread)(A,v);CHKERRQ(ierr);
  CHKMEMQ;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseRestoreArrayRead"
/*@C
   MatDenseRestoreArrayRead - Restores the matrix after MatDenseGetArrayRead() has been called.

   Not Collective

   Input Parameter:
+  A - the matrix
-  v - the location of the values

   Level: advanced

.seealso: MatDenseGetArrayRead()
@*/

PetscErrorCode  MatDenseRestoreArrayRead(MatDense A,const PetscScalar *v[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidPointer(v,2);
#if defined(PETSC_USE_DEBUG)
  CHKMEMQ;
#endif
  if (!A->ops->restorearrayread) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDense type %s",((PetscObject)A)->type_name);
  ierr = (*A->ops->restorearrayread)(A,v);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "MatDenseDuplicate"
/*@
   MatDenseDuplicate - Duplicates a matrix.

   Collective on MatDense

   Input Parameters:
+  A - the matrix
-  op - MATDENSE_COPY_VALUES copies the numerical values in the matrix, MATDENSE_DO_NOT_COPY_VALUES avoids that, and MATDENSE_MAKE_SIBLING makes A and M share the same data.

   Output Parameter:
.  M - pointer to place new matrix

   Level: intermediate

   Concepts: matrices^duplicating

.seealso: MatDenseCopy()
@*/

PetscErrorCode  MatDenseDuplicate(MatDense A,MatDenseDuplicateOption op,MatDense *M)
{
  PetscErrorCode ierr;
  MatDense       B;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidPointer(M,3);
  ierr = PetscCheckMatDenseForRead(A);CHKERRQ(ierr);

  *M = 0;
  if (!A->ops->duplicate) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not written for this matrix type");
  ierr = PetscLogEventBegin(MATDENSE_Duplicate,A,0,0,0);CHKERRQ(ierr);
  ierr = (*A->ops->duplicate)(A,op,M);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(MATDENSE_Duplicate,A,0,0,0);CHKERRQ(ierr);
  B = *M;
  if (op == MATDENSE_MAKE_SIBLING) {
    B->ld = A->ld;
    B->is_allocated = A->is_allocated;
    B->is_hermitian = A->is_hermitian;
    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);
}


#undef __FUNCT__  
#define __FUNCT__ "MatDenseAreSiblings"
/*@C
   MatDenseAreSiblings - return whether one matrix comes from MatDenseDuplicte
   of the other with the MATDENSE_MAKE_SIBLING option.

   Not Collective

   Input Parameter:
.  A,B - the matrices

   Output Parameter:
.  aresiblings - whether the matrices are siblings

   Level: advanced

.seealso: MatDenseDuplicate()
@*/

PetscErrorCode  MatDenseAreSiblings(MatDense A,MatDense B,PetscBool *aresiblings)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidType(A,1);
  PetscValidHeaderSpecific(B,MATDENSE_CLASSID,2);
  PetscValidType(B,2);
  PetscValidPointer(aresiblings,3);
  if (A->ops->aresiblings != B->ops->aresiblings)
    *aresiblings = PETSC_FALSE;
  else if (A->ops->aresiblings) {
    ierr = (*A->ops->aresiblings)(A,B,aresiblings);CHKERRQ(ierr);
  } else {
    *aresiblings = A->data == B->data ? PETSC_TRUE : PETSC_FALSE;
  }
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "MatDenseSetFromOptions"
/*@
   MatDenseSetFromOptions - Creates a matrix where the type is determined
   from the options database. The default matrix type is
   VEC, using the routines MatDenseCreateVec() if
   you do not select a type in the options database.

   Collective on MatDense

   Input Parameter:
.  A - the matrix

   Options Database Keys:
+    -mat_type vec     - VEC type, uses MatDenseCreateVec()
-    -mat_type magma   - MAGMA type, uses MatDenseCreateMAGMA()

   Even More Options Database Keys:
   See the manpages for particular formats (e.g., MatDenseCreateVec())
   for additional format-specific options.

   Level: beginner

.keywords: matrix, create

.seealso: MatDenseCreateVec(), MatDenseCreateMAGMA()
@*/

PetscErrorCode  MatDenseSetFromOptions(MatDense A)
{
  PetscErrorCode ierr;
  const char     *deft = MATDENSEBASIC;
  char           type[256];
  PetscBool      flg,opb;
  PetscInt       opi;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);

  ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"Dense matrix options","MatDense");CHKERRQ(ierr);
    ierr = PetscOptionsList("-matdense_type","Dense matrix type","MatDenseSetType",MatDenseList,deft,type,256,&flg);CHKERRQ(ierr);
    if (flg) {
      ierr = MatDenseSetType(A,type);CHKERRQ(ierr);
    } else if (!((PetscObject)A)->type_name) {
      ierr = MatDenseSetType(A,deft);CHKERRQ(ierr);
    }
    ierr = PetscOptionsInt("-matdense_buffersize","buffersize(KB) for matrix-matrix multiplication","",A->matmult_buffersize/1024,&opi,&flg);CHKERRQ(ierr);
    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);
    }

    /* process any options handlers added with PetscObjectAddOptionsHandler() */
    ierr = PetscObjectProcessOptionsHandlers((PetscObject)A);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseSetUpPreallocation"
/*@
   MatDenseSetUpPreallocation

   Collective on MatDense

   Input Parameter:
.  B - the matrix

   Level: beginner

.keywords: matrix, create
@*/

PetscErrorCode  MatDenseSetUpPreallocation(MatDense B)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(B,MATDENSE_CLASSID,1);
  if (B->is_allocated == PETSC_FALSE && B->ops->setuppreallocation) {
    ierr = PetscInfo(B,"Warning not preallocating matrix storage\n");CHKERRQ(ierr);
    ierr = PetscLogEventBegin(MATDENSE_SetUpPreallocation,B,0,0,0);CHKERRQ(ierr);
    ierr = (*B->ops->setuppreallocation)(B);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(MATDENSE_SetUpPreallocation,B,0,0,0);CHKERRQ(ierr);
    B->is_allocated = PETSC_TRUE;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MatDenseSetType"
/*@C
   MatDenseSetType - Builds matrix object for a particular matrix type

   Collective on MatDense

   Input Parameters:
+  mat      - the matrix object
-  matype   - matrix type

   Options Database Key:
.  -matdense_type  <method> - Sets the type; use -help for a list
    of available methods (for instance, seqaij)

  Level: intermediate

.keywords: MatDense, MatDenseType, set, method

.seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat
@*/

PetscErrorCode  MatDenseSetType(MatDense mat, const MatDenseType matype)
{
  PetscErrorCode ierr,(*r)(MatDense);
  PetscBool      sametype;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MATDENSE_CLASSID,1);

  ierr = PetscTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr);
  if (sametype) PetscFunctionReturn(0);

  ierr =  PetscFListFind(MatDenseList,((PetscObject)mat)->comm,matype,PETSC_TRUE,(void(**)(void))&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatDense type given: %s",matype);
 
  /* free the old data structure if it existed */
  if (mat->ops->destroy) {
    ierr = (*mat->ops->destroy)(mat);CHKERRQ(ierr);
    mat->ops->destroy = PETSC_NULL;
  }
  mat->is_allocated = PETSC_FALSE;

  /* create the new data structure */
  ierr = (*r)(mat);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode MatDenseView_Default_ASCII(MatDense xin,PetscViewer viewer);

#undef __FUNCT__  
#define __FUNCT__ "MatDenseView"
/*@C
   MatDenseView - Visualizes a dense matrix object.

   Collective on MatDense

   Input Parameters:
+  mat - the matrix
-  viewer - visualization context

  Notes:
  The available visualization contexts include
+    PETSC_VIEWER_STDOUT_SELF - standard output (default)
.    PETSC_VIEWER_STDOUT_WORLD - synchronized standard
        output where only the first processor opens
        the file.  All other processors send their
        data to the first processor to print.
-     PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure

   The user can open alternative visualization contexts with
+    PetscViewerASCIIOpen() - Outputs matrix to a specified file
.    PetscViewerBinaryOpen() - Outputs matrix in binary to a
         specified file; corresponding input uses MatLoad()
.    PetscViewerDrawOpen() - Outputs nonzero matrix structure to
         an X window display
-    PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
         Currently only the sequential dense and AIJ
         matrix types support the Socket viewer.

   The user can call PetscViewerSetFormat() to specify the output
   format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
   PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
+    PETSC_VIEWER_DEFAULT - default, prints matrix contents
.    PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
.    PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
.    PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
         format common among all matrix types
.    PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
         format (which is in many cases the same as the default)
.    PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
         size and structure (not the matrix entries)
.    PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
         the matrix structure

   Level: beginner

   Concepts: matrices^viewing

.seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
          PetscViewerSocketOpen(), PetscViewerBinaryOpen()
@*/

PetscErrorCode  MatDenseView(MatDense mat,PetscViewer viewer)
{
  PetscErrorCode    ierr;
  PetscInt          rows,cols;
  PetscBool         iascii;
  PetscViewerFormat format;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MATDENSE_CLASSID,1);
  PetscValidType(mat,1);
  if (!viewer) {
    ierr = PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);CHKERRQ(ierr);
  }
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
  PetscCheckSameComm(mat,1,viewer,2);
  ierr = PetscCheckMatDenseForRead(mat);CHKERRQ(ierr);

  ierr = PetscLogEventBegin(MATDENSE_View,mat,viewer,0,0);CHKERRQ(ierr);
  ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
  if (iascii) {
    ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);  
    if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
      ierr = PetscObjectPrintClassNamePrefixType((PetscObject)mat,viewer,"Matrix Object");CHKERRQ(ierr);
      ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
      ierr = MatDenseGetSizes(mat,PETSC_NULL,PETSC_NULL,&rows,&cols);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(viewer,"rows=%D, cols=%D\n",rows,cols);CHKERRQ(ierr);
    }
  }
  if (mat->ops->view) {
    ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
    ierr = (*mat->ops->view)(mat,viewer);CHKERRQ(ierr);
    ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
  } 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) {
    ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);  
    if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
      ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
    }
  }
  ierr = PetscLogEventEnd(MATDENSE_View,mat,viewer,0,0);CHKERRQ(ierr);
  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
  MatDenseFinalizePackage - This function destroys everything in the Slepc interface to the MatDense package. It is
  called from SlepcFinalize().

  Level: developer

.seealso: SlepcFinalize()
@*/

PetscErrorCode MatDenseFinalizePackage(void)
{
  PetscFunctionBegin;
  MatDensePackageInitialized = PETSC_FALSE;
  MatDenseList               = 0;
  MatDenseRegisterAllCalled  = PETSC_FALSE;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatDenseInitializePackage"
/*@C
  MatDenseInitializePackage - This function initializes everything in the EPSDense package. It is called
  from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to MatDenseCreate()
  when using static libraries.

  Input Parameter:
  path - The dynamic library path, or PETSC_NULL

  Level: developer

.seealso: SlepcInitialize()
@*/

PetscErrorCode MatDenseInitializePackage(const char *path) {
  char           logList[256];
  char           *className;
  PetscBool      opt;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (MatDensePackageInitialized) PetscFunctionReturn(0);
  MatDensePackageInitialized = PETSC_TRUE;
  /* Register Classes */
  ierr = PetscClassIdRegister("MatDense",&MATDENSE_CLASSID);CHKERRQ(ierr);
  /* Register Constructors */
  ierr = MatDenseRegisterAll(path);CHKERRQ(ierr);
  /* Register Events */
  ierr = PetscLogEventRegister("MatDenseDuplicate",MATDENSE_CLASSID,&MATDENSE_Duplicate);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDensePrealloc",MATDENSE_CLASSID,&MATDENSE_SetUpPreallocation);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseMatMult",MATDENSE_CLASSID,&MATDENSE_MatMult);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseBlasMatMult",MATDENSE_CLASSID,&MATDENSE_BlasMatMult);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseCopy",MATDENSE_CLASSID,&MATDENSE_Copy);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseAXPY",MATDENSE_CLASSID,&MATDENSE_AXPY);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseView",MATDENSE_CLASSID,&MATDENSE_View);CHKERRQ(ierr);
  ierr = PetscLogEventRegister("MatDenseConvert",MATDENSE_CLASSID,&MATDENSE_Convert);CHKERRQ(ierr);
  /* Process info exclusions */
  ierr = PetscOptionsGetString(PETSC_NULL,"-info_exclude",logList,256,&opt);CHKERRQ(ierr);
  if (opt) {
    ierr = PetscStrstr(logList,"matdense",&className);CHKERRQ(ierr);
    if (className) {
      ierr = PetscInfoDeactivateClass(MATDENSE_CLASSID);CHKERRQ(ierr);
    }
  }
  /* Process summary exclusions */
  ierr = PetscOptionsGetString(PETSC_NULL,"-log_summary_exclude",logList,256,&opt);CHKERRQ(ierr);
  if (opt) {
    ierr = PetscStrstr(logList,"matdense",&className);CHKERRQ(ierr);
    if (className) {
      ierr = PetscLogEventDeactivateClass(MATDENSE_CLASSID);CHKERRQ(ierr);
    }
  }
  ierr = PetscRegisterFinalize(MatDenseFinalizePackage);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatDenseGetType"
/*@C
   MatDenseGetType - Gets the MatDense type as a string from the MatDense object.

   Not Collective

   Input Parameter:
.  A - the dense matrix

   Output Parameter:
.  name - name of MatDense method

   Level: intermediate

.seealso: MatDenseSetType()
@*/

PetscErrorCode MatDenseGetType(MatDense A,const MatDenseType *type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MATDENSE_CLASSID,1);
  PetscValidPointer(type,2);
  *type = ((PetscObject)A)->type_name;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatDenseRegister"
/*@C
  MatDenseRegister - See MatDenseRegisterDynamic()

  Level: advanced
@*/

PetscErrorCode MatDenseRegister(const char *sname,const char *path,const char *name,PetscErrorCode (*function)(MatDense))
{
  PetscErrorCode ierr;
  char           fullname[PETSC_MAX_PATH_LEN];

  PetscFunctionBegin;
  ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr);
  ierr = PetscFListAdd(&MatDenseList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatDenseRegisterDestroy"
/*@
   MatDenseRegisterDestroy - Frees the list of MatDense methods that were
   registered by MatDenseRegisterDynamic().

   Not Collective

   Level: advanced

.seealso: MatDenseRegisterDynamic(), MatDenseRegisterAll()
@*/

PetscErrorCode MatDenseRegisterDestroy(void)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscFListDestroy(&MatDenseList);CHKERRQ(ierr);
  MatDenseRegisterAllCalled = PETSC_FALSE;
  PetscFunctionReturn(0);
}