I don't fully understand what the empty KSP is doing either, because it seems to me that if a KSP object was properly destroyed, KSPCreate() should have to be called again before it is used again, right? So it must not be fully destroyed or else the second linear system solve would fail, but something is happening to it such that the symbolic factorization is redone.
Like I said, the fix I implemented seems to work (even if it doesn't seem necessary), so I don't know that it really needs any more attention. However, I was able to reproduce the behavior in the example you sent me. If you're curious, running the program as normal should behave the same way as the original example. Running with the option "-destroy_eps" destroys and recreates the EPS object, and running with the additional option "-empty_ksp" gives the EPS object an empty KSP before it is destroyed. At least when I ran it on my machine, I got 1 symbolic factorization for the original program, 2 with the "-destroy_eps" option, and 1 with both "-destroy_eps" and "-empty_ksp". Thanks, Darin ________________________________________ From: Jose E. Roman [[email protected]] Sent: Wednesday, November 02, 2016 3:29 AM To: Peetz, Darin T Cc: [email protected] Subject: Re: [petsc-users] Provide Matrix Factorization to EPS for Generalized Eigenvalue Problem > El 1 nov 2016, a las 20:42, Peetz, Darin T <[email protected]> escribió: > > Thank you for providing the example. I have managed to get it working > properly; however, after looking at the example, I have one more question. > > First, let me provide a little background on how the code is set up. I have > one function that takes care of solving the linear system of equations, and a > separate one that takes care of the eigenvalue problem. The KSP object for > the linear system is passed between both functions, and kept for the duration > of the program to save some of the cost of factorizing the matrix each time. > I have been keeping the EPS object local to the function that takes care of > the eigenvalue problem since it is only relevant within that function, and I > haven't seen a need to reuse the same EPS object in every iteration (I'm > assuming the cost of creating the EPS object itself is relatively small). > The reason for the extra symbolic factorizations seems to have been that when > I destroyed the EPS object before it goes out of scope, it also destroys the > underlying KSP object that I provided. > > Basically, I'm wondering if my assumption about the cost of creating the EPS > object is correct. Is it okay to recreate the EPS object every time the > eigenvalue function is called, or does this create a significant overhead? > It seems that if I try to reuse the same EPS object I have to make sure I'm > reusing the same A-matrix as well, or else the ST object redoes both the > symbolic and numeric factorizations (so the same K matrix is numerically > factorized twice in each iteration) each time. Thus I would be able to save > a significant amount of memory if I could let the EPS object and A-matrix > exist in a local scope. > > As a side note, I worked around the KSP destruction issue by providing an > empty (uninitialized) KSP object to the EPS object after calling EPSSolve and > extracting the eigenvalues/eigenvectors. This is functional but seems like a > bad/crude approach. Is there a more elegant way to "free" the KSP object > from the EPS object once it is no longer needed, or should I stick with this > approach? > > Thanks a lot for your help, > Darin STSetKSP() increases the reference count of the KSP object, so when EPS is destroyed, the KSP object remains. I don't understand why you need to provide an empty KSP. Regarding your question about recreating EPS, yes the cost can be considered small. Jose
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SLEPc - Scalable Library for Eigenvalue Problem Computations Copyright (c) 2002-2016, 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/>. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ static char help[] = "Generalized Symmetric eigenproblem.\n\n" "The problem is Ax = lambda Bx, with:\n" " A = Laplacian operator in 2-D\n" " B = diagonal matrix\n\n" "The command line options are:\n" " -n <n>, where <n> = number of grid subdivisions in x dimension.\n" " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n"; #include <slepceps.h> #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc,char **argv) { Mat A,B; /* matrices */ EPS eps; /* eigenproblem solver context */ ST st; /* spectral transformation context */ Vec x,b; KSP ksp; PC pc; PetscInt N,n=10,m,Istart,Iend,II,i,j; PetscReal norm; PetscBool flag,terse; PetscErrorCode ierr; PetscBool destroy, empty; ierr = SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);CHKERRQ(ierr); if (!flag) m=n; N = n*m; ierr = PetscOptionsHasName(NULL,NULL,"-terse",&terse);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized Symmetric Eigenproblem, N=%D (%Dx%D grid)\n\n",N,n,m);CHKERRQ(ierr); /// Flags dealing with optional destruction of EPS object PetscOptionsHasName(NULL,NULL,"-destroy_eps",&destroy); PetscOptionsHasName(NULL,NULL,"-empty_ksp",&empty); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute the matrices that define the eigensystem, Ax=kBx - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); ierr = MatSetUp(A);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); ierr = MatSetFromOptions(B);CHKERRQ(ierr); ierr = MatSetUp(B);CHKERRQ(ierr); ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); for (II=Istart;II<Iend;II++) { i = II/n; j = II-i*n; if (i>0) { ierr = MatSetValue(A,II,II-n,-1.0,INSERT_VALUES);CHKERRQ(ierr); } if (i<m-1) { ierr = MatSetValue(A,II,II+n,-1.0,INSERT_VALUES);CHKERRQ(ierr); } if (j>0) { ierr = MatSetValue(A,II,II-1,-1.0,INSERT_VALUES);CHKERRQ(ierr); } if (j<n-1) { ierr = MatSetValue(A,II,II+1,-1.0,INSERT_VALUES);CHKERRQ(ierr); } ierr = MatSetValue(A,II,II,4.0,INSERT_VALUES);CHKERRQ(ierr); ierr = MatSetValue(B,II,II,4.0,INSERT_VALUES);CHKERRQ(ierr); } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatCreateVecs(B,&x,&b);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create a KSP with coefficient matrix B - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve first linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = VecSetRandom(b,NULL);CHKERRQ(ierr); ierr = KSPSetOperators(ksp,B,B);CHKERRQ(ierr); ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of x= %g\n",(double)norm);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create the eigensolver - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr); ierr = EPSSetOperators(eps,A,B);CHKERRQ(ierr); ierr = EPSSetProblemType(eps,EPS_GHEP);CHKERRQ(ierr); ierr = EPSGetST(eps,&st);CHKERRQ(ierr); ierr = STSetKSP(st,ksp);CHKERRQ(ierr); ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve first eigensystem - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = EPSSolve(eps);CHKERRQ(ierr); if (terse) { ierr = EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);CHKERRQ(ierr); } else { ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); ierr = EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Optionally destroy EPS object - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (empty) { KSP empty; ierr = EPSGetST(eps,&st);CHKERRQ(ierr); ierr = STSetKSP(st,empty);CHKERRQ(ierr); } if (destroy || empty) { ierr = EPSDestroy(&eps);CHKERRQ(ierr); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Modify B matrix - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ for (II=Istart;II<Iend;II++) { i = II/n; j = II-i*n; ierr = MatSetValue(B,II,II,2.0+i,INSERT_VALUES);CHKERRQ(ierr); } ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve second linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = VecSetRandom(b,NULL);CHKERRQ(ierr); ierr = KSPSetOperators(ksp,B,B);CHKERRQ(ierr); ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of x= %g\n",(double)norm);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Optionally recreate the eigensolver - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (destroy || empty) { ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr); } ierr = EPSSetOperators(eps,A,B);CHKERRQ(ierr); if (destroy || empty) { ierr = EPSSetProblemType(eps,EPS_GHEP);CHKERRQ(ierr); ierr = EPSGetST(eps,&st);CHKERRQ(ierr); ierr = STSetKSP(st,ksp);CHKERRQ(ierr); ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve second eigensystem - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = EPSSolve(eps);CHKERRQ(ierr); if (terse) { ierr = EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);CHKERRQ(ierr); } else { ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); ierr = EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); } ierr = EPSDestroy(&eps);CHKERRQ(ierr); ierr = KSPDestroy(&ksp);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = MatDestroy(&B);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = SlepcFinalize(); return ierr; }
