Subversion Repositories slepc-dev

/*

   Basic implementation of 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 <slepcblaslapack.h>

typedef struct {

  PetscScalar *v;

  PetscBool   user_allocated;

  PetscInt    siblings;

} MatDense_Basic;

static PetscErrorCode  MatDenseCreateBasic_Private(MatDense A, MatDense_Basic *b,PetscScalar *data);

#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_,*vC;

  const PetscScalar

                 *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);

  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;

  mB = Bt == PETSC_TRUE ? B->n : B->m;

  nB = Bt == PETSC_TRUE ? B->m : B->n;

  ldA = A->ld;

  ldB = B->ld;

  ldC = C->ld;

  qA = At == PETSC_TRUE ? T : N;

  qB = Bt == PETSC_TRUE ? T : N;

  if (mA == 1 && nA == 1 && nB == 1) {

    if (!At && !Bt)     *vC =           *vA  *           *vB;

    else if (At && !Bt) *vC = PetscConj(*vA) *           *vB;

    else if (!At && Bt) *vC =           *vA  * PetscConj(*vB);

    else if (At && Bt)  *vC = PetscConj(*vA) * PetscConj(*vB);

    ierr = PetscLogFlops(1);CHKERRQ(ierr);

  /* 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 !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 !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 !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); }

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

  }

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

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseBlasCopy_Basic"

PetscErrorCode  MatDenseBlasCopy_Basic(PetscInt n,MatDense A,PetscInt dA,PetscInt iA,MatDense B,PetscInt dB,PetscInt iB)

{

  PetscErrorCode ierr;

  PetscScalar    *vB,*vB_;

  const PetscScalar

                 *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);

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

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseBlasAXPY_Basic"

PetscErrorCode  MatDenseBlasAXPY_Basic(PetscInt n,PetscScalar a,MatDense X,PetscInt dX,PetscInt iX,MatDense Y,PetscInt dY,PetscInt iY)

{

  PetscErrorCode ierr;

  PetscScalar    *vY,*vY_;

  const PetscScalar

                 *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);

  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 = PetscLogFlops(2*n);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseSetExplicit_Basic"

PetscErrorCode  MatDenseSetExplicit_Basic(MatDense A)

{

  PetscErrorCode ierr;

  PetscScalar    *vA,*vA_;

  PetscInt       i;

  PetscBLASInt   n,incX,one=1;

  PetscFunctionBegin;

  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;

      BLAScopy_(&n,&vA[i],&incX,&vA[A->ld*i+i+1],&one);

    }

  } else if (MatDenseIsImplicitTriangular(A)) {

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

      ierr = PetscMemzero(&vA[A->ld*i+i+1],(A->n-i-1)*sizeof(PetscScalar));CHKERRQ(ierr);

    }

  }

  A->is_impl = PETSC_FALSE;

  ierr = MatDenseRestoreArray(A,&vA_);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__

#define __FUNCT__ "MatDenseDestroy_Basic"

PetscErrorCode  MatDenseDestroy_Basic(MatDense A)

{

  PetscErrorCode ierr;

  MatDense_Basic *data = (MatDense_Basic*)A->data;

  PetscFunctionBegin;

  if (--data->siblings <= 0) {

    if (!data->user_allocated) { ierr = PetscFree(data->v);CHKERRQ(ierr); }

    ierr = PetscFree(data);CHKERRQ(ierr);

  }

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseBlasAXPY_*_*_C","",PETSC_NULL);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseBlasCopy_*_*_C","",PETSC_NULL);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseMatMatMult_*_*_C","",PETSC_NULL);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseGetArray_Basic"

PetscErrorCode  MatDenseGetArray_Basic(MatDense A,PetscScalar *v[])

{

  MatDense_Basic *data = (MatDense_Basic*)A->data;

  PetscFunctionBegin;

  *v = data->v;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseRestoreArray_Basic"

PetscErrorCode  MatDenseRestoreArray_Basic(MatDense A,PetscScalar *v[])

{

  PetscFunctionBegin;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseGetArrayRead_Basic"

PetscErrorCode  MatDenseGetArrayRead_Basic(MatDense A,const PetscScalar *v[])

{

  MatDense_Basic *data = (MatDense_Basic*)A->data;

  PetscFunctionBegin;

  *v = data->v;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseRestoreArrayRead_Basic"

PetscErrorCode  MatDenseRestoreArrayRead_Basic(MatDense A,const PetscScalar *v[])

{

  PetscFunctionBegin;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseDuplicate_Basic"

PetscErrorCode  MatDenseDuplicate_Basic(MatDense A,MatDenseDuplicateOption op,MatDense *M)

{

  MatDense_Basic *data = (MatDense_Basic*)A->data;

  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = MatDenseCreate(((PetscObject)A)->comm,M);CHKERRQ(ierr);

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

  if (op == MATDENSE_COPY_VALUES) {

    ierr = MatDenseSetUpPreallocation(*M);CHKERRQ(ierr);

    ierr = MatDenseCopy(A,*M);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseSetFromOptions_Basic"

PetscErrorCode  MatDenseSetFromOptions_Basic(MatDense A)

{

  PetscFunctionBegin;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "MatDenseSetUpPreallocation_Basic"

PetscErrorCode  MatDenseSetUpPreallocation_Basic(MatDense B)

{

  PetscErrorCode ierr;

  MatDense_Basic *data = (MatDense_Basic*)B->data;

  PetscFunctionBegin;

  if (data->v) PetscFunctionReturn(0);

  ierr = PetscMalloc(B->ld*B->Nmax*sizeof(PetscScalar),&data->v);CHKERRQ(ierr);

  ierr = PetscLogObjectMemory(B,B->ld*B->Nmax*sizeof(PetscScalar));CHKERRQ(ierr);

  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,

       MatDenseRestoreArray_Basic,

       MatDenseGetArrayRead_Basic,

       MatDenseRestoreArrayRead_Basic,

       MatDenseDuplicate_Basic,

  /* 5 */

       0/*MatDenseAXPY_Basic*/,

       0/*MatDenseCopy_Basic*/,

       MatDenseSetExplicit_Basic,

       MatDenseBlasAXPY_Basic,

       MatDenseBlasCopy_Basic,

  /* 10 */

       MatDenseDestroy_Basic,

       0/*MatDenseView_Basic*/,

       MatDenseSetFromOptions_Basic,

       0/*MatDenseMatMult_Basic*/,

       MatDenseBlasMatMult_Basic,

  /* 15 */

       0/*MatDenseSetSizes_Basic*/,

       MatDenseSetUpPreallocation_Basic,

       0/*MatDenseAreSiblings_Basic*/,

       MatDenseSerialize_Basic,

       MatDenseDeserialize_Basic,

};

#undef __FUNCT__  

#define __FUNCT__ "MatDenseCreateBasic_Private"

static PetscErrorCode  MatDenseCreateBasic_Private(MatDense A, MatDense_Basic *b, PetscScalar *v)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  if (b) { b->siblings++; }

  else {

    ierr = PetscNewLog(A,MatDense_Basic,&b);CHKERRQ(ierr);

    b->v = v;

    b->user_allocated = v?PETSC_TRUE:PETSC_FALSE;

    b->siblings = 1;

  }

  A->data = (void*)b;

  ierr = PetscMemcpy(A->ops,&MatDenseOps_Basic,sizeof(struct _MatDenseOps));CHKERRQ(ierr);

  ierr = PetscObjectChangeTypeName((PetscObject)A,MATDENSEBASIC);CHKERRQ(ierr);

  A->is_allocated = b->v?PETSC_TRUE:PETSC_FALSE;

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseBlasAXPY_*_*_C","MatDenseBlasAXPY_Basic",MatDenseBlasAXPY_Basic);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseBlasCopy_*_*_C","MatDenseBlasCopy_Basic",MatDenseBlasCopy_Basic);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatDenseBlasMatMult_*_*_C","MatDenseMatMult_Basic",MatDenseBlasMatMult_Basic);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

/*MC

   MATDENSEBASIC - A basic implementation of dense matrix using standard BLAS

                   and LAPACK functions.

  Level: advanced

.seealso: MatDenseCreateBasic()

M*/

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "MatDenseCreate_Basic"

PetscErrorCode  MatDenseCreate_Basic(MatDense A)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = MatDenseCreateBasic_Private(A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

EXTERN_C_END

#undef __FUNCT__  

#define __FUNCT__ "MatDenseCreateBasic"

/*@C

   MatDenseCreateBasic - Creates a new dense matrix that uses standard BLAS

                         and LAPACK functions.

  Collective on MPI_Comm

   Input Parameters:

+  comm - MPI communicator

.  M - maximum number of rows

.  N - maximum number of columns

.  m - current number of rows

.  n - current number of columns

-  data - optional location of the matrix data. Set to PETSC_NULL for PETSc

   to control its memory allocation.

   Output Parameter:

.  A - the matrix

   Level: advanced

.keywords: matrix, mat, create

.seealso: MATDENSEBASIC

@*/

PetscErrorCode  MatDenseCreateBasic(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscScalar *data,MatDense *A)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = MatDenseCreate(comm,A);CHKERRQ(ierr);

  ierr = MatDenseSetMaxSizes(*A,M,N);CHKERRQ(ierr);

  ierr = MatDenseCreateBasic_Private(*A,PETSC_NULL,data);CHKERRQ(ierr);

  ierr = MatDenseSetSizes(*A,0,0,m,n);CHKERRQ(ierr);

  ierr = MatDenseSetUpPreallocation(*A);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}