Subversion Repositories slepc-dev

/*

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