Subversion Repositories slepc-dev

/*

   Basic PS routines

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2011, Universitat 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 <slepc-private/psimpl.h>      /*I "slepcps.h" I*/

PetscFList       PSList = 0;

PetscBool        PSRegisterAllCalled = PETSC_FALSE;

PetscClassId     PS_CLASSID = 0;

PetscLogEvent    PS_Solve = 0,PS_Sort = 0,PS_Vectors = 0,PS_Other = 0;

static PetscBool PSPackageInitialized = PETSC_FALSE;

const char       *PSMatName[PS_NUM_MAT] = {"A","B","C","T","Q","X","Y","U","VT","W"};

#undef __FUNCT__  

#define __FUNCT__ "PSFinalizePackage"

/*@C

   PSFinalizePackage - This function destroys everything in the Slepc interface

   to the PS package. It is called from SlepcFinalize().

   Level: developer

.seealso: SlepcFinalize()

@*/

PetscErrorCode PSFinalizePackage(void)

{

  PetscFunctionBegin;

  PSPackageInitialized = PETSC_FALSE;

  PSList               = 0;

  PSRegisterAllCalled  = PETSC_FALSE;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSInitializePackage"

/*@C

  PSInitializePackage - This function initializes everything in the PS package.

  It is called from PetscDLLibraryRegister() when using dynamic libraries, and

  on the first call to PSCreate() when using static libraries.

  Input Parameter:

  path - The dynamic library path, or PETSC_NULL

  Level: developer

.seealso: SlepcInitialize()

@*/

PetscErrorCode PSInitializePackage(const char *path)

{

  char             logList[256];

  char             *className;

  PetscBool        opt;

  PetscErrorCode   ierr;

  PetscFunctionBegin;

  if (PSPackageInitialized) PetscFunctionReturn(0);

  PSPackageInitialized = PETSC_TRUE;

  /* Register Classes */

  ierr = PetscClassIdRegister("Projected system",&PS_CLASSID);CHKERRQ(ierr);

  /* Register Constructors */

  ierr = PSRegisterAll(path);CHKERRQ(ierr);

  /* Register Events */

  ierr = PetscLogEventRegister("PSSolve",PS_CLASSID,&PS_Solve);CHKERRQ(ierr);

  ierr = PetscLogEventRegister("PSSort",PS_CLASSID,&PS_Sort);CHKERRQ(ierr);

  ierr = PetscLogEventRegister("PSVectors",PS_CLASSID,&PS_Vectors);CHKERRQ(ierr);

  ierr = PetscLogEventRegister("PSOther",PS_CLASSID,&PS_Other);CHKERRQ(ierr);

  /* Process info exclusions */

  ierr = PetscOptionsGetString(PETSC_NULL,"-info_exclude",logList,256,&opt);CHKERRQ(ierr);

  if (opt) {

    ierr = PetscStrstr(logList,"ps",&className);CHKERRQ(ierr);

    if (className) {

      ierr = PetscInfoDeactivateClass(PS_CLASSID);CHKERRQ(ierr);

    }

  }

  /* Process summary exclusions */

  ierr = PetscOptionsGetString(PETSC_NULL,"-log_summary_exclude",logList,256,&opt);CHKERRQ(ierr);

  if (opt) {

    ierr = PetscStrstr(logList,"ps",&className);CHKERRQ(ierr);

    if (className) {

      ierr = PetscLogEventDeactivateClass(PS_CLASSID);CHKERRQ(ierr);

    }

  }

  ierr = PetscRegisterFinalize(PSFinalizePackage);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSCreate"

/*@C

   PSCreate - Creates a PS context.

   Collective on MPI_Comm

   Input Parameter:

.  comm - MPI communicator

   Output Parameter:

.  newps - location to put the PS context

   Level: beginner

   Note:

   PS objects are not intended for normal users but only for

   advanced user that for instance implement their own solvers.

.seealso: PSDestroy(), PS

@*/

PetscErrorCode PSCreate(MPI_Comm comm,PS *newps)

{

  PS             ps;

  PetscInt       i;

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidPointer(newps,2);

  ierr = PetscHeaderCreate(ps,_p_PS,struct _PSOps,PS_CLASSID,-1,"PS","Projected System","PS",comm,PSDestroy,PSView);CHKERRQ(ierr);

  *newps      = ps;

  ps->state   = PS_STATE_RAW;

  ps->method  = 0;

  ps->nmeth   = 1;

  ps->compact = PETSC_FALSE;

  ps->ld      = 0;

  ps->l       = 0;

  ps->n       = 0;

  ps->k       = 0;

  for (i=0;i<PS_NUM_MAT;i++) {

    ps->mat[i]  = PETSC_NULL;

    ps->rmat[i] = PETSC_NULL;

  }

  ps->work   = PETSC_NULL;

  ps->rwork  = PETSC_NULL;

  ps->iwork  = PETSC_NULL;

  ps->lwork  = 0;

  ps->lrwork = 0;

  ps->liwork = 0;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSetOptionsPrefix"

/*@C

   PSSetOptionsPrefix - Sets the prefix used for searching for all

   PS options in the database.

   Logically Collective on PS

   Input Parameters:

+  ps - the projected system context

-  prefix - the prefix string to prepend to all PS option requests

   Notes:

   A hyphen (-) must NOT be given at the beginning of the prefix name.

   The first character of all runtime options is AUTOMATICALLY the

   hyphen.

   Level: advanced

.seealso: PSAppendOptionsPrefix()

@*/

PetscErrorCode PSSetOptionsPrefix(PS ps,const char *prefix)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  ierr = PetscObjectSetOptionsPrefix((PetscObject)ps,prefix);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSAppendOptionsPrefix"

/*@C

   PSAppendOptionsPrefix - Appends to the prefix used for searching for all

   PS options in the database.

   Logically Collective on PS

   Input Parameters:

+  ps - the projected system context

-  prefix - the prefix string to prepend to all PS option requests

   Notes:

   A hyphen (-) must NOT be given at the beginning of the prefix name.

   The first character of all runtime options is AUTOMATICALLY the hyphen.

   Level: advanced

.seealso: PSSetOptionsPrefix()

@*/

PetscErrorCode PSAppendOptionsPrefix(PS ps,const char *prefix)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  ierr = PetscObjectAppendOptionsPrefix((PetscObject)ps,prefix);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__

#define __FUNCT__ "PSGetOptionsPrefix"

/*@C

   PSGetOptionsPrefix - Gets the prefix used for searching for all

   PS options in the database.

   Not Collective

   Input Parameters:

.  ps - the projected system context

   Output Parameters:

.  prefix - pointer to the prefix string used is returned

   Notes: On the fortran side, the user should pass in a string 'prefix' of

   sufficient length to hold the prefix.

   Level: advanced

.seealso: PSSetOptionsPrefix(), PSAppendOptionsPrefix()

@*/

PetscErrorCode PSGetOptionsPrefix(PS ps,const char *prefix[])

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidPointer(prefix,2);

  ierr = PetscObjectGetOptionsPrefix((PetscObject)ps,prefix);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSetType"

/*@C

   PSSetType - Selects the type for the PS object.

   Logically Collective on PS

   Input Parameter:

+  ps   - the projected system context

-  type - a known type

   Level: advanced

.seealso: PSGetType()

@*/

PetscErrorCode PSSetType(PS ps,const PSType type)

{

  PetscErrorCode ierr,(*r)(PS);

  PetscBool      match;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidCharPointer(type,2);

  ierr = PetscTypeCompare((PetscObject)ps,type,&match);CHKERRQ(ierr);

  if (match) PetscFunctionReturn(0);

  ierr =  PetscFListFind(PSList,((PetscObject)ps)->comm,type,PETSC_TRUE,(void (**)(void))&r);CHKERRQ(ierr);

  if (!r) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PS type %s",type);

  ierr = PetscMemzero(ps->ops,sizeof(struct _PSOps));CHKERRQ(ierr);

  ierr = PetscObjectChangeTypeName((PetscObject)ps,type);CHKERRQ(ierr);

  ierr = (*r)(ps);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSGetType"

/*@C

   PSGetType - Gets the PS type name (as a string) from the PS context.

   Not Collective

   Input Parameter:

.  ps - the projected system context

   Output Parameter:

.  name - name of the projected system

   Level: advanced

.seealso: PSSetType()

@*/

PetscErrorCode PSGetType(PS ps,const PSType *type)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidPointer(type,2);

  *type = ((PetscObject)ps)->type_name;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSetMethod"

/*@

   PSSetMethod - Selects the method to be used to solve the problem.

   Logically Collective on PS

   Input Parameter:

+  ps   - the projected system context

-  meth - an index indentifying the method

   Level: advanced

.seealso: PSGetMethod()

@*/

PetscErrorCode PSSetMethod(PS ps,PetscInt meth)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidLogicalCollectiveInt(ps,meth,2);

  if (meth<0) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The method must be a non-negative integer");

  ps->method = meth;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSGetMethod"

/*@

   PSGetMethod - Gets the method currently used in the PS.

   Not Collective

   Input Parameter:

.  ps - the projected system context

   Output Parameter:

.  meth - identifier of the method

   Level: advanced

.seealso: PSSetMethod()

@*/

PetscErrorCode PSGetMethod(PS ps,PetscInt *meth)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidPointer(meth,2);

  *meth = ps->method;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSetCompact"

/*@

   PSSetCompact - Switch to compact storage of matrices.

   Logically Collective on PS

   Input Parameter:

+  ps   - the projected system context

-  comp - a boolean flag

   Notes:

   Compact storage is used in some PS types such as PSHEP when the matrix

   is tridiagonal. This flag can be used to indicate whether the user

   provides the matrix entries via the compact form (the tridiagonal PS_MAT_T)

   or the non-compact one (PS_MAT_A).

   The default is PETSC_FALSE.

   Level: advanced

.seealso: PSGetCompact()

@*/

PetscErrorCode PSSetCompact(PS ps,PetscBool comp)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidLogicalCollectiveBool(ps,comp,2);

  ps->compact = comp;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSGetCompact"

/*@

   PSGetCompact - Gets the compact storage flag.

   Not Collective

   Input Parameter:

.  ps - the projected system context

   Output Parameter:

.  comp - the flag

   Level: advanced

.seealso: PSSetCompact()

@*/

PetscErrorCode PSGetCompact(PS ps,PetscBool *comp)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidPointer(comp,2);

  *comp = ps->compact;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSetFromOptions"

/*@

   PSSetFromOptions - Sets PS options from the options database.

   Collective on PS

   Input Parameters:

.  ps - the projected system context

   Notes:  

   To see all options, run your program with the -help option.

   Level: beginner

@*/

PetscErrorCode PSSetFromOptions(PS ps)

{

  PetscErrorCode ierr;

  PetscInt       meth;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  if (!PSRegisterAllCalled) { ierr = PSRegisterAll(PETSC_NULL);CHKERRQ(ierr); }

  /* Set default type (we do not allow changing it with -ps_type) */

  if (!((PetscObject)ps)->type_name) {

    ierr = PSSetType(ps,PSNHEP);CHKERRQ(ierr);

  }

  ierr = PetscOptionsBegin(((PetscObject)ps)->comm,((PetscObject)ps)->prefix,"Projecte System (PS) Options","PS");CHKERRQ(ierr);

    meth = 0;

    ierr = PetscOptionsInt("-ps_method","Method to be used for the projected system","PSSetMethod",ps->method,&meth,PETSC_NULL);CHKERRQ(ierr);

    ierr = PSSetMethod(ps,meth);CHKERRQ(ierr);

    ierr = PetscObjectProcessOptionsHandlers((PetscObject)ps);CHKERRQ(ierr);

  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSView"

/*@C

   PSView - Prints the PS data structure.

   Collective on PS

   Input Parameters:

+  ps - the projected system context

-  viewer - optional visualization context

   Note:

   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.

   The user can open an alternative visualization context with

   PetscViewerASCIIOpen() - output to a specified file.

   Level: beginner

.seealso: PetscViewerASCIIOpen()

@*/

PetscErrorCode PSView(PS ps,PetscViewer viewer)

{

  PetscBool         isascii;

  const char        *state;

  PetscViewerFormat format;

  PetscErrorCode    ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  if (!viewer) viewer = PETSC_VIEWER_STDOUT_(((PetscObject)ps)->comm);

  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);

  PetscCheckSameComm(ps,1,viewer,2);

  ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);

  if (isascii) {

    ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);

    ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ps,viewer,"PS Object");CHKERRQ(ierr);

    if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {

      switch (ps->state) {

        case PS_STATE_RAW:          state = "raw"; break;

        case PS_STATE_INTERMEDIATE: state = "intermediate"; break;

        case PS_STATE_CONDENSED:    state = "condensed"; break;

        case PS_STATE_SORTED:       state = "sorted"; break;

        default: SETERRQ(((PetscObject)ps)->comm,1,"Wrong value of ps->state");

      }

      ierr = PetscViewerASCIIPrintf(viewer,"  current state: %s\n",state);CHKERRQ(ierr);

      ierr = PetscViewerASCIIPrintf(viewer,"  dimensions: ld=%d, n=%d, l=%d, k=%d\n",ps->ld,ps->n,ps->l,ps->k);CHKERRQ(ierr);

    }

    if (ps->ops->view) {

      ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);

      ierr = (*ps->ops->view)(ps,viewer);CHKERRQ(ierr);

      ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);

    }

  } else SETERRQ1(((PetscObject)ps)->comm,1,"Viewer type %s not supported for PS",((PetscObject)viewer)->type_name);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSAllocate"

/*@

   PSAllocate - Allocates memory for internal storage or matrices in PS.

   Logically Collective on PS

   Input Parameters:

+  ps - the projected system context

-  ld - leading dimension (maximum allowed dimension for the matrices)

   Level: advanced

.seealso: PSGetLeadingDimension(), PSSetDimensions()

@*/

PetscErrorCode PSAllocate(PS ps,PetscInt ld)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ps,PS_CLASSID,1);

  PetscValidLogicalCollectiveInt(ps,ld,2);

  if (ld<1) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Leading dimension should be at least one");

  ps->ld = ld;

  ierr = (*ps->ops->allocate)(ps,ld);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSAllocateMat_Private"

PetscErrorCode PSAllocateMat_Private(PS ps,PSMatType m)

{

  PetscInt       sz;

  PetscErrorCode ierr;

  PetscFunctionBegin;

  if (m==PS_MAT_T) sz = 3*ps->ld*sizeof(PetscScalar);

  else sz = ps->ld*ps->ld*sizeof(PetscScalar);

  if (ps->mat[m]) { ierr = PetscFree(ps->mat[m]);CHKERRQ(ierr); }

  else { ierr = PetscLogObjectMemory(ps,sz);CHKERRQ(ierr); }

  ierr = PetscMalloc(sz,&ps->mat[m]);CHKERRQ(ierr);

  ierr = PetscMemzero(ps->mat[m],sz);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSAllocateMatReal_Private"

PetscErrorCode PSAllocateMatReal_Private(PS ps,PSMatType m)

{

  PetscInt       sz;

  PetscErrorCode ierr;

  PetscFunctionBegin;

  if (m==PS_MAT_T) sz = 3*ps->ld*sizeof(PetscReal);

  else sz = ps->ld*ps->ld*sizeof(PetscReal);

  if (ps->rmat[m]) { ierr = PetscFree(ps->rmat[m]);CHKERRQ(ierr); }

  else { ierr = PetscLogObjectMemory(ps,sz);CHKERRQ(ierr); }

  ierr = PetscMalloc(sz,&ps->rmat[m]);CHKERRQ(ierr);

  ierr = PetscMemzero(ps->rmat[m],sz);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSAllocateWork_Private"

PetscErrorCode PSAllocateWork_Private(PS ps,PetscInt s,PetscInt r,PetscInt i)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  if (s>ps->lwork) {

    ierr = PetscFree(ps->work);CHKERRQ(ierr);

    ierr = PetscMalloc(s*sizeof(PetscScalar),&ps->work);CHKERRQ(ierr);

    ierr = PetscLogObjectMemory(ps,(s-ps->lwork)*sizeof(PetscScalar));CHKERRQ(ierr);

    ps->lwork = s;

  }

  if (r>ps->lrwork) {

    ierr = PetscFree(ps->rwork);CHKERRQ(ierr);

    ierr = PetscMalloc(r*sizeof(PetscReal),&ps->rwork);CHKERRQ(ierr);

    ierr = PetscLogObjectMemory(ps,(r-ps->lrwork)*sizeof(PetscReal));CHKERRQ(ierr);

    ps->lrwork = r;

  }

  if (i>ps->liwork) {

    ierr = PetscFree(ps->iwork);CHKERRQ(ierr);

    ierr = PetscMalloc(i*sizeof(PetscBLASInt),&ps->iwork);CHKERRQ(ierr);

    ierr = PetscLogObjectMemory(ps,(i-ps->liwork)*sizeof(PetscBLASInt));CHKERRQ(ierr);

    ps->liwork = i;

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSViewMat_Private"

PetscErrorCode PSViewMat_Private(PS ps,PetscViewer viewer,PSMatType m)

{

  PetscErrorCode    ierr;

  PetscInt          i,j;

  PetscScalar       *v;

  PetscViewerFormat format;

#if defined(PETSC_USE_COMPLEX)

  PetscBool         allreal = PETSC_TRUE;

#endif

  PetscFunctionBegin;

  ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);

  if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0);

  ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);

#if defined(PETSC_USE_COMPLEX)

  /* determine if matrix has all real values */

  v = ps->mat[m];

  for (i=0;i<ps->n;i++)

    for (j=0;j<ps->n;j++)

      if (PetscImaginaryPart(v[i+j*ps->ld])) { allreal = PETSC_FALSE; break; }

#endif

  if (format == PETSC_VIEWER_ASCII_MATLAB) {

    ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",ps->n,ps->n);CHKERRQ(ierr);

    ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",PSMatName[m]);CHKERRQ(ierr);

  } else {

    ierr = PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",PSMatName[m]);CHKERRQ(ierr);

  }

  for (i=0;i<ps->n;i++) {

    v = ps->mat[m]+i;

    for (j=0;j<ps->n;j++) {

#if defined(PETSC_USE_COMPLEX)

      if (allreal) {

        ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr);

      } else {

        ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr);

      }

#else

      ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr);

#endif

      v += ps->ld;

    }

    ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);

  }

  if (format == PETSC_VIEWER_ASCII_MATLAB) {

    ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr);

  }

  ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);

  ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSortEigenvaluesReal_Private"

PetscErrorCode PSSortEigenvaluesReal_Private(PS ps,PetscInt n,PetscReal *eig,PetscInt *perm,PetscErrorCode (*comp_func)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void *comp_ctx)

{

  PetscErrorCode ierr;

  PetscScalar    re;

  PetscInt       i,j,result,tmp;

  PetscFunctionBegin;

  for (i=0;i<n;i++) perm[i] = i;

  /* insertion sort */

  for (i=1;i<n;i++) {

    re = eig[perm[i]];

    j = i-1;

    ierr = (*comp_func)(re,0.0,eig[perm[j]],0.0,&result,comp_ctx);CHKERRQ(ierr);

    while (result<=0 && j>=0) {

      tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;

      if (j>=0) {

        ierr = (*comp_func)(re,0.0,eig[perm[j]],0.0,&result,comp_ctx);CHKERRQ(ierr);

      }

    }

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSPermuteColumns_Private"

PetscErrorCode PSPermuteColumns_Private(PS ps,PetscInt n,PSMatType m,PetscInt *perm)

{

  PetscInt    i,j,k,p,ld;

  PetscScalar *Q,rtmp;

  PetscFunctionBegin;

  ld = ps->ld;

  Q  = ps->mat[PS_MAT_Q];

  for (i=0;i<n;i++) {

    p = perm[i];

    if (p != i) {

      j = i + 1;

      while (perm[j] != i) j++;

      perm[j] = p; perm[i] = i;

      /* swap eigenvectors i and j */

      for (k=0;k<n;k++) {

        rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;

      }

    }

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSReset"

/*@

   PSReset - Resets the PS context to the initial state.

   Collective on PS

   Input Parameter:

.  ps - the projected system context