Subversion Repositories slepc-dev

/*

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

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

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

*/

static char help[] = "Tests matdense interface.\n\n";

#include <slepcmatdense.h>

#undef __FUNCT__

#define __FUNCT__ "main"

int main( int argc, char **argv )

{

  MatDense       A,B,C;        /* matrices */

  PetscScalar    *v;

  PetscInt       n=45,m,i,j,its,M,reps=1;

  PetscBool      flag;

  Vec            *vecs0,*vecs1,t,*vecso;

  PetscLogStage  saxpy,smult,supdate,supdateo;

  PetscErrorCode ierr;

  SlepcInitialize(&argc,&argv,(char*)0,help);

  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);

  ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);CHKERRQ(ierr);

  if(!flag) m=n;

  M = PetscMax(m,n);

  ierr = PetscOptionsGetInt(PETSC_NULL,"-r",&reps,&flag);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"\nDense matrices of size %Dx%D\n\n",m,n);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     Initialize vectors

     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = VecCreate(PETSC_COMM_WORLD,&t);CHKERRQ(ierr);

  ierr = VecSetSizes(t,PETSC_DECIDE,m);CHKERRQ(ierr);

  ierr = VecSetFromOptions(t);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(t,n,&vecso);CHKERRQ(ierr);

  ierr = SlepcVecSetTemplate(t);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(t,n,&vecs0);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(t,n,&vecs1);CHKERRQ(ierr);

  ierr = VecGetLocalSize(t,&m);CHKERRQ(ierr);

  ierr = VecDestroy(&t);CHKERRQ(ierr);

  for (i=0;i<n;i++) { ierr = VecSetRandom(vecs0[i],PETSC_NULL);CHKERRQ(ierr); }

  for (i=0;i<n;i++) { ierr = VecSetRandom(vecs1[i],PETSC_NULL);CHKERRQ(ierr); }

  for (i=0;i<n;i++) { ierr = VecSetRandom(vecso[i],PETSC_NULL);CHKERRQ(ierr); }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     Create matrices

     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = MatDenseCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);

  ierr = MatDenseSetMaxSizes(A,m,n);CHKERRQ(ierr);

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

  ierr = MatDenseSetFromOptions(A);CHKERRQ(ierr);

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

  ierr = MatDenseDuplicate(A,MATDENSE_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr);

  ierr = MatDenseSetFromOptions(B);CHKERRQ(ierr);

  ierr = MatDenseSetUpPreallocation(B);CHKERRQ(ierr);

  ierr = MatDenseCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);

  ierr = MatDenseSetMaxSizes(C,M,M);CHKERRQ(ierr);

  ierr = MatDenseSetFromOptions(C);CHKERRQ(ierr);

  ierr = MatDenseSetUpPreallocation(C);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     Initialize matrices

     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);

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

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

      v[m*j+i] = 1.0;

    }

  }

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

  ierr = MatDenseGetArray(B,&v);CHKERRQ(ierr);

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

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

      v[m*j+i] = i & 1 ? 1.0 : -1.0;

    }

  }

  ierr = MatDenseRestoreArray(B,&v);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     Test operations

     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = PetscLogStageRegister("AXPY",&saxpy);CHKERRQ(ierr);

  ierr = PetscLogStagePush(saxpy);CHKERRQ(ierr);

  for (i=0; i<reps; i++) { ierr = MatDenseAXPY(B,1.0,A);CHKERRQ(ierr); }

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

    for (j=0; j<n; j++) { ierr = VecAXPY(vecs0[j],1.0,vecs1[j]);CHKERRQ(ierr); }

  }

  PetscLogStagePop();

  ierr = PetscLogStageRegister("MatMult",&smult);CHKERRQ(ierr);

  ierr = PetscLogStagePush(smult);CHKERRQ(ierr);

  ierr = MatDenseSetSizes(C,0,0,n,n);CHKERRQ(ierr);

  for (i=0; i<reps; i++) { ierr = MatDenseMatMult(C,1.0,0.0,A,PETSC_TRUE,B,PETSC_FALSE);CHKERRQ(ierr); }

  ierr = MatDenseGetArray(C,&v);CHKERRQ(ierr);

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

    for (j=0; j<n; j++) { ierr = VecMDot(vecs0[j],n,vecs1,&v[n*j]);CHKERRQ(ierr); }

  }

  ierr = MatDenseRestoreArray(C,&v);CHKERRQ(ierr);

  PetscLogStagePop();

  ierr = PetscLogStageRegister("Update",&supdate);CHKERRQ(ierr);

  ierr = PetscLogStagePush(supdate);CHKERRQ(ierr);

  ierr = MatDenseSetSizes(C,0,0,n,n);CHKERRQ(ierr);

  ierr = MatDenseGetArray(C,&v);CHKERRQ(ierr);

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

    ierr = SlepcUpdateVectors(n,vecs0,0,n,v,m,PETSC_FALSE);CHKERRQ(ierr);

    for (j=0; j<n; j++) { ierr = VecAXPY(vecs0[j],1.0,vecs1[j]);CHKERRQ(ierr); }

  }

  ierr = MatDenseRestoreArray(C,&v);CHKERRQ(ierr);

  PetscLogStagePop();

  ierr = PetscLogStageRegister("Update Orig.",&supdateo);CHKERRQ(ierr);

  ierr = PetscLogStagePush(supdateo);CHKERRQ(ierr);

  ierr = MatDenseSetSizes(C,0,0,n,n);CHKERRQ(ierr);

  ierr = MatDenseGetArray(C,&v);CHKERRQ(ierr);

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

    ierr = SlepcUpdateVectors(n,vecso,0,n,v,m,PETSC_FALSE);CHKERRQ(ierr);

    for (j=0; j<n; j++) { ierr = VecAXPY(vecso[j],1.0,vecs1[j]);CHKERRQ(ierr); }

  }

  ierr = MatDenseRestoreArray(C,&v);CHKERRQ(ierr);

  PetscLogStagePop();

  ierr = MatDenseDestroy(&A);CHKERRQ(ierr);

  ierr = MatDenseDestroy(&B);CHKERRQ(ierr);

  ierr = MatDenseDestroy(&C);CHKERRQ(ierr);

  ierr = VecDestroyVecs(n,&vecs0);CHKERRQ(ierr);

  ierr = VecDestroyVecs(n,&vecs1);CHKERRQ(ierr);

  ierr = VecDestroyVecs(n,&vecso);CHKERRQ(ierr);

  ierr = SlepcFinalize();CHKERRQ(ierr);

  return 0;

}