hello
i am trying to solve generalized eigenvalue problem using eps solver. i
have a monte carlo loop in the program which has to run for about 1000
times. matrix sizes are 15 X 15 as a dummy problem and 23040 X23040 in the
real problem.
as i run the program i see that eps takes increasingly more amount of time
for computation with respect to iterations and number of processors too.
i have used all the options to the best of my knowledge and need your help
to improve the program in any way possible.
i am attaching .c file of mycode with this email.
I would be extremely greatful to you for your valuable suggestions.
*Siddhesh M Godbole*
5th year Dual Degree,
Civil Eng & Applied Mech.
IIT Madras
static char help[]="eps";
#include <petscmat.h>
#include <petscvec.h>
#include <../src/mat/impls/sbaij/seq/sbaij.h>
#include <petscblaslapack.h>
#include <petscksp.h>
#include <math.h>
#include <petsc-private/vecimpl.h>
#include <petsc-private/taoimpl.h>
#include "../src/tao/pde_constrained/impls/lcl/lcl.h"
#include <slepceps.h>
#include "mpi.h"
extern PetscErrorCode CkEigenSolutions(PetscInt*,Mat*,PetscReal*,Vec*,PetscInt*,PetscInt*,PetscReal*);
extern PetscErrorCode FormElementStiffness_idx0(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscScalar*);
extern PetscErrorCode FormElementStiffness(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscScalar*);
extern PetscErrorCode FormElementMass_idx0(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);
extern PetscErrorCode FormElementMass(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);
extern PetscErrorCode FormRandomElementStiffness(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
Mat K,M,DM,DK,ME,KE,Phi,TPhi,P,disp,vel,acc,f,invM,Kcap,a,b,invKcap,U,mmtop,kktop,a2DM;
Vec Iv,*Cv,xr,xi,delpe;
Vec random,P0,P1,P2,Pans,dKcap,delp,delu,delv,dela,v1,v1a3,acc1a5,v2,*av1,bacc1,v1a4,acc1a6,a1delu,a2delu,disp1,disp2,acc1,acc2/*(intial condition vector of acceleration)*/;
char file1[PETSC_MAX_PATH_LEN],file2[PETSC_MAX_PATH_LEN]; /* input file name */
PetscBool flg,flgB=PETSC_FALSE;
PetscErrorCode ierr;
PetscRandom rnd;
PetscBool preload=PETSC_TRUE,isSymmetric,evecs,isHermitian,isGeneralized,isPositive,*aK;
PetscScalar *arrayK,*arrayM,*evecs_array,kr,ki,re,im,*axr;
PetscMPIInt size,rank;
PetscInt i,j,z,k,nev,il,iu,r,trial=1,change;
PetscInt maxit,its,lits,nconv,nini=1,ncon=0,S=100000,Dim=100000;
PetscLogStage stages[2];
PetscReal vl,vu,abstol=1.e-8;
PetscScalar Matlab_E,freq[Dim],total[Dim],average[Dim],REV[trial],Xr[Dim];
PetscInt t,/*n=5,m=2*/id[1],idx[6],IDX[3],index[1],idxmcsn[Dim],idxmcsm[1],idxrandom[trial],idx1[Dim],idxnt[1000],rows[Dim],one[1],istart,iend,nrows;
/*S = (n+1)*(m);*/
ST st;
EPS eps;
EPSType type;
PetscScalar Ke[16], Me[16],K0[4],M0[4];
PetscScalar width = 25.65/1000; /*25.65 mm */
PetscScalar depth = 3.25/1000; /* 3.25 mm */
PetscScalar Area = width*depth, rho = 2662; /*2662 kg/m3*/
PetscScalar Tlel =.31; /*length of the beam*/
PetscScalar lel = Tlel/5; /*length of each element*/
PetscScalar I = width*depth*depth*depth/12; /*moment of inertia*/
PetscScalar E = 50000000000,En; /*Modulus of elasticity*/
PetscViewer fd1,fd2,fd3,viewer;
PetscReal tol;
/*parameters used in Newmarks method*/
PetscScalar value[1000],*aP0,*aDM,*aDK,*aa,*ab,*aKcap,*aav1,*abacc1,*aP2,*adelp,*adelu,*av1a3,*av1a4,*aacc1a5,*aacc1a6,*adelv,*adela,*aa1delu,*aa2delu,*adisp1,*adisp2,*arv1,*arv2,*aacc1,*aacc2;
PetscScalar ainvKcap[S],*adKcap;
PetscScalar gaama=1/2,beta=1/4,dt=.01,tt=1000,a1=200,a2= 40000,a3=400,a4= 2,a5= 2,a6=0;
SlepcInitialize(&argc,&args,(char*)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
/*Setting up socket connection for Petsc-Matlab interfacing*/
/*fd1 = PETSC_VIEWER_SOCKET_WORLD;
/*fd2 = PETSC_VIEWER_SOCKET_2;
/*fd3 = PETSC_VIEWER_SOCKET_3;
PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd1);*/
/*PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd3);
/*if matrices are to be read from the file*/
/*ierr = PetscOptionsGetString(NULL,"-file1",file1,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading Real 'Mass' matrix from a binary file...\n");CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file1,FILE_MODE_READ,&viewer);CHKERRQ(ierr);*/
/*ierr = MatCreate(PETSC_COMM_WORLD,&mmtop);CHKERRQ(ierr);
ierr = MatLoad(mmtop,fd1);CHKERRQ(ierr);
/*ierr = MatView(mmtop,fd1);CHKERRQ(ierr);
/*ierr = MatSetFromOptions(Mfile);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);*/
/*ierr = PetscOptionsGetString(NULL,"-file2",file2,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading Real 'Stiffness' matrix from a binary file...\n");CHKERRQ(ierr);
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file2,FILE_MODE_READ,&fd);CHKERRQ(ierr);*/
/*PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd2); */
/*ierr = MatCreate(PETSC_COMM_WORLD,&kktop);CHKERRQ(ierr);
ierr = MatLoad(kktop,fd1);CHKERRQ(ierr);
/*ierr = MatView(kktop,fd1);CHKERRQ(ierr);
/*ierr = MatSetFromOptions(Kfile);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);*/
PetscRandomCreate(PETSC_COMM_WORLD,&rnd);
PetscRandomSetType(rnd,PETSCRAND48);
PetscRandomSetFromOptions(rnd);
for(i=0;i<trial;i++) { PetscRandomGetValue(rnd,&Matlab_E); REV[i]=Matlab_E-0.5; }
/* for(i=0;i<trial/5;i++) printf("from processor %d: %f %f %f %f %f \n ",rank,REV[i],REV[trial/5+i],REV[2*trial/5+i],REV[3*trial/5+i],REV[4*trial/5+i]);*/
PetscRandomDestroy(&rnd);
for(i=0;i<S;i++) {freq[i]=0;total[i]=0;average[i]=0;}
/*_____________________________________________________________________ Creating mass matrix_______________________________________________________*/
ierr = MatCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);
ierr = MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,S,S);CHKERRQ(ierr);
ierr = MatSetFromOptions(M);CHKERRQ(ierr);
ierr = MatSetUp(M);CHKERRQ(ierr);
/*
Assemble matrix
*/
for (i=0; i<5; i++)
{
/* node numbers for the four corners of element */
if(i==0)
{ ierr = FormElementMass_idx0(rho,lel,width,depth,M0);
IDX[0]=0;
IDX[1]=1;
IDX[2]=2;
ierr = MatSetValues(M,3,IDX,3,IDX,M0,ADD_VALUES);CHKERRQ(ierr);}
else { ierr = FormElementMass(rho,lel,width,depth,Me);
idx[0] = 3*(i-1);
idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1; idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
ierr = MatSetValues(M,6,idx,6,idx,Me,ADD_VALUES);CHKERRQ(ierr);}
}
ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*MatView(M,fd1);
/* printf("Mass Matrix is" "/n" "/n");
*/
/*PetscLogStagePush(stages[0]);*/
change=0;
/* ******************************************************************LOOP FOR MONTE CARLO SIMULATION STARTS****************************************************** */
for(z=0;z<trial;z++)
{
/*printf( " \n \n %d th trial on process %d out of %d processes is working fine! \n", z,rank+1, size);
/* _________________________________________________________________________________Stiffness matrices ___________________________________________________________________*/
ierr = MatCreate(PETSC_COMM_WORLD,&K);CHKERRQ(ierr);
ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,S,S);CHKERRQ(ierr);
ierr = MatSetFromOptions(K);CHKERRQ(ierr);
ierr = MatSetUp(K);CHKERRQ(ierr);
/*
Assemble matrix
*/
for (i=0; i<5; i++)
{
if(i==0)
{
ierr = FormElementStiffness_idx0(E,Area,lel,I,K0);
IDX[0]=2*i;
IDX[1]=2*i+1;
IDX[2]=2*i+2;
ierr = MatSetValues(K,3,IDX,3,IDX,K0,ADD_VALUES);CHKERRQ(ierr);
}
else if(i<4)
{
ierr = FormElementStiffness(E,Area,lel,I,Ke);CHKERRQ(ierr);
idx[0] = 3*(i-1);
idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1; idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
ierr = MatSetValues(K,6,idx,6,idx,Ke,ADD_VALUES);CHKERRQ(ierr);}
else
{
En = E*(1 + REV[z]/10);
ierr = FormRandomElementStiffness(En,Area,lel,I,Ke);
idx[0] = 3*(i-1);
idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1; idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
ierr = MatSetValues(K,6,idx,6,idx,Ke,ADD_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* printf("Mass Matrix is" "/n" "/n");
MatView(K,viewer);*/
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and set various options
/*
Create eigensolver context
*/
ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
/*
Set operators. In this case, it is a generalized eigenvalue problem
*/
ierr = EPSSetOperators(eps,K,M);CHKERRQ(ierr); /*==========================Do necessary changes with Kfile and Mfile according to the location from which they are taken======*/
MatGetVecs(K,NULL,&xr);
/* VecSet(xr,100);*/
/*MatGetVecs(K,NULL,&xi); */ /*=========Do necessary changes=======*/
/*
If the user provided initial guesses or constraints, pass them here
*/
/*VecLoad(Iv,fd1);*/
ierr = EPSSetInitialSpace(eps,1,&xr);CHKERRQ(ierr);
/* ierr = EPSSetDeflationSpace(eps,ncon,Cv);CHKERRQ(ierr);*/
/*
Set solver parameters at runtime
*/
ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the eigensystem - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*set problem as generalized hermitian*/
EPSSetProblemType(eps,EPS_GHEP);
ierr = EPSIsGeneralized(eps,&isGeneralized);CHKERRQ(ierr);
ierr = EPSIsHermitian(eps,&isHermitian);CHKERRQ(ierr);
ierr = EPSIsPositive(eps,&isPositive);CHKERRQ(ierr);
ierr = EPSSolve(eps);CHKERRQ(ierr);
/*ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr);*/
ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
/*printf(" converged eigenpairs are %d \n",nconv);*/
/* setting up the eigenvector matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&Phi);CHKERRQ(ierr);
ierr = MatSetSizes(Phi,PETSC_DECIDE,PETSC_DECIDE,Dim,nconv);CHKERRQ(ierr);
ierr = MatSetFromOptions(Phi);CHKERRQ(ierr);
ierr = MatSetUp(Phi);CHKERRQ(ierr);
/*Save eigenvectors, if requested in Phi matrix*/
/*initializing Phi matrix*/
for (i=0;i<nconv;i++)
{
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
/*ierr = PetscPrintf(PETSC_COMM_WORLD,"eigenvalues are: %12f \n",sqrt((double)kr)/(2*3.1415926535));CHKERRQ(ierr);*/
VecGetArray(xr,&axr);
for(j=0;j<Dim;j++){idx1[j]=j;}
MatSetValues(Phi,Dim,&idx1,1,&i,axr,ADD_VALUES);
VecRestoreArray(xr,&axr);
/*ierr = padmatrix(xr,S,nconv,Phi,i);CHKERRQ(ierr);*/
re = kr;
im = ki;
total[i]= total[i] + kr;
}
/*ierr = VecDestroy(&xr);CHKERRQ(ierr);*/
ierr = MatAssemblyBegin(Phi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Phi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*ierr = MatView(Phi,fd1);CHKERRQ(ierr);
/* ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);*/
ierr = PetscFree(xi);CHKERRQ(ierr);
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);
/*_________________________________________________________________________________________________________________________________
Solving system dynamics
_________________________________________________________________________________________________________________________________*/
ierr = MatTranspose(Phi,MAT_INITIAL_MATRIX,&TPhi);CHKERRQ(ierr); /* B = A^T */
ierr = MatMatMult(M,Phi,MAT_INITIAL_MATRIX,1,&ME);CHKERRQ(ierr); /*ME = M*Phi*/
ierr = MatMatMult(TPhi,ME,MAT_INITIAL_MATRIX,1,&DM);CHKERRQ(ierr); /*DM = TPhi*ME*/
ierr = MatMatMult(K,Phi,MAT_INITIAL_MATRIX,1,&KE);CHKERRQ(ierr); /*ME = K*Phi*/
ierr = MatMatMult(TPhi,KE,MAT_INITIAL_MATRIX,1,&DK);CHKERRQ(ierr); /*DM = TPhi*KE*/
/*ierr = PetscPrintf(PETSC_COMM_WORLD,"Decoupled mass matrix is\n");CHKERRQ(ierr);*/
/* MatView(DM,fd1);
/*ierr = PetscPrintf(PETSC_COMM_WORLD,"Decoupled stiffness matrix is\n");CHKERRQ(ierr);
MatView(DK,PETSC_VIEWER_STDOUT_WORLD);*/
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
/*ierr = MatDestroy(&Phi);CHKERRQ(ierr);*/
ierr = MatDestroy(&KE);CHKERRQ(ierr);
ierr = MatDestroy(&ME);CHKERRQ(ierr);
/* ************************************************************************** LOOP FOR MONTE CARLO ENDS******************************************************/
}
for(i=0;i<nconv;i++) {
average[i]=total[i]/trial;
PetscPrintf(PETSC_COMM_WORLD,"average eigenvalues are :%f \n",sqrt(average[i])/(2*3.1415926535));
}
/* Free work space. */
ierr = MatDestroy(&M);CHKERRQ(ierr);
/*ierr = PetscFinalize();*/
ierr = SlepcFinalize();
return 0;
}
