El 28/05/2014, a las 19:27, Giacomo Mulas escribió: > Hello. > > After some time stuck doing less exciting stuff, I got back to trying to use > slepc for science. I am trying to use the relatively new functionality for > arbitrary selection of the eigenpairs to be found, and I would want some > support to understand if I am doing things correctly. I apologise with Jose > Roman for disappearing when I should have helped testing this, it was not my > choice. But now I am back at it. > > In particular, I want to obtain the eigenvectors with the maximum projection > in a given subspace (defined by a number of normalised vectors, let's call > them targets). > > I don't know in advance how many eigenpairs must be determined: I want to > obtain enough eigenvectors that the projection of all targets in the space > of these eigenvectors is very nearly identical to the targets themselves. > > So my strategy, so far, is the following: > > 1) create and set up the matrix H to be diagonalised > > ierr = MatCreate(mixpars->slepc_comm, &H); CHKERRQ(ierr); > ierr = MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, statesinlist, > statesinlist); > CHKERRQ(ierr); > ierr = MatSetFromOptions(H);CHKERRQ(ierr); > ierr = MatSetUp(H);CHKERRQ(ierr); > ... > ... some MatSetValue(H,...); > ... > ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > ierr = MatGetVecs(H, &xr, PETSC_NULL);CHKERRQ(ierr); > ierr = MatGetVecs(H, &xi, PETSC_NULL);CHKERRQ(ierr); > > 2) create the eps > > ierr = EPSCreate(mixpars->slepc_comm,&eps);CHKERRQ(ierr); > ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr); > ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); > ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE); > ierr = EPSSetFromOptions(eps); CHKERRQ(ierr); > ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr); > ierr = EPSSetFromOptions(eps); CHKERRQ(ierr); > ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr); > > 3) set the eps to use an arbitrary selection function > > /* every time I solve, I want to find one eigenvector, and it must be > the one with the largest component along the target state */ > ierr = EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE); > ierr = EPSSetArbitrarySelection(eps, computeprojection(), > (void *) &targetcompindex); > ierr = EPSSetDimensions(eps, 1, PETSC_IGNORE, PETSC_IGNORE); > CHKERRQ(ierr); > > 4) run a loop over the target vectors, and iteratively call EPSSolve until > each target vector is completely contained in the space of the > eigenvectors found. Before each call to EPSSolve, I set the initial guess > equal to the target vector, and set the deflation space to be the set of > eigenvectors found so far. After each call to EPSSolve, I add the new > eigenvectors to the deflation space one by one, and check if the target > state is (nearly) fully contained in the eigenvectors space. If yes, I > move on to the next target state and so on. > > 5) free everything, destroy eps, matrices, vectors etc. > > > I have some questions about the above: > > 1) should it work in principle, or am I getting it all wrong?
I don't see much problem. > > 2) should I destroy and recreate the eps after each call to EPSSolve and > before next call? Or, since the underlying matrix is always the same, can I > just call EPSSetInitialSpace(), EPSSetDeflationSpace(), update the internal > parameter to be passed to the arbitrary selection function and I can call > again EPSSolve? No need to recreate the solver. The only thing is EPSSetDeflationSpace() - I would suggest calling EPSRemoveDeflationSpace() and then EPSSetDeflationSpace() again with the extended set of vectors. Do not call EPSSetDeflationSpace() with a single vector every time. > > 3) Since what I want is going on to find one eigenpair at a time of the same > problem until some condition is fulfilled, is there a way in which I can > achieve this without setting it up again and again every time? Can I specify > an arbitrary function that is called by EPSSolve to decide whether enough > eigenpairs were computed or not, instead of doing it in this somewhat > awkward manner? No. > > 4) more technical: since I add vectors one by one to the deflation space, to > begin with I allocate a Vec *Cv with PetscMalloc(statesinlist*sizeof(Vec *), > &Cv); > where statesinlist is the size of the problem, hence the maximum > hypothetical size of the deflation space. I would prefer to allocate this > dynamically, enlarging it as needed. Is there something like realloc() in > PETSC/SLEPC? > > Thanks in advance, bye > Giacomo Mulas > > -- > _________________________________________________________________ > > Giacomo Mulas <[email protected]> > _________________________________________________________________ > > INAF - Osservatorio Astronomico di Cagliari > via della scienza 5 - 09047 Selargius (CA) > > tel. +39 070 71180244 > mob. : +39 329 6603810 > _________________________________________________________________ > > "When the storms are raging around you, stay right where you are" > (Freddy Mercury) > _________________________________________________________________
