Another beginner's querry:
I worked out an example that "combines"
mat/examples/tests/ex145.c
/ksp/ksp/examples/tutorials/ex2.c
/ksp/ksp/examples/tutorials/ex30.c
Please see attached code sample. The code compiles but does not solve the
system correctly. I would like to know what changes I would need to
implement to get it working right. Thanks!
P.S. Feel free to include the (corrected) code sample in the test/tutorial
folders of the development version.
On Sun, Apr 12, 2015 at 7:47 PM, Hong <[email protected]> wrote:
> You can either use elemental directly, or use petsc-elemental interface.
> An example can be found at
> ~petsc/src/mat/examples/tests/ex145.c
>
> You may use petsc KSP interface instead. I just modified
> ~petsc/src/ksp/ksp/examples/tutorials/ex2.c
> so this example can be run with elemental with runtime options
> mpiexec -n 3 ./ex2 -pc_type lu -pc_factor_mat_solver_package elemental
> -mat_type elemental
> Norm of error 2.81086e-15 iterations 1
>
> Please using petsc-dev (master branch) for petsc-elemental interface.
>
> Hong
>
> On Sun, Apr 12, 2015 at 6:57 PM, Preyas Shah <[email protected]>
> wrote:
>
>> Hi,
>>
>> I have been recently investigating the use of Petsc for solving a PDE
>> related to my research and web search suggests that I should use petsc with
>> elemental.
>>
>> So far, I was required to solve a matrix equation Ax=b where A was dense
>> (with number of non zeros =0) but the size of the matrix was order
>> 5000x5000. I employed the standard serial LU solver from Gnu Scientific
>> Library and obtained a decent runtime that served my needs.
>>
>> Now I am investigating the same problem in a particularly singular limit
>> of one parameters in my PDE. As a result, to obtain grid convergence on the
>> physical domain, I am forced to go to sizes of A beyond 30000x30000. I am
>> trying to find a good library that can solve such dense systems in
>> **parallel**. Petsc says its capable of doing dense linear algebra but my
>> web search hasn't shown me any examples where dense equations are solved in
>> **parallel**. A webpage showing a minimum working example would be enough.
>> Or any other advice :)
>>
>> Thanks for your time!
>> ~Preyas
>>
>>
>
--
~
Preyas
Doctoral Researcher, MechE, Stanford Univ.
*"Can a man still be brave if he's afraid?".**"That is the only time he can
be brave."-- George R R Martin*
*"To learn who rules over you, simply find out who you are not allowed to
criticize." -- Voltaire*
*"First they ignore you, then they laugh at you, then they fight you, then
you win" -- Gandhi*
#include <petscksp.h>
#include <iostream>
#include "mpi.h"
#define RUNKSP
#ifdef RUNKSP
int main(int argc, char* argv[]){
/*-----------------------------------------------------------------------------
* Declare vars and objects
*-----------------------------------------------------------------------------*/
Vec rhs,sol,ex_sol;
Mat MM;
KSP ksp;
PC pc;
PetscReal norm;
PetscInt i,j,nrows,ncols,its,tb=1000,tt;
tt = tb;
PetscErrorCode ierr;
const PetscInt *rows;
const PetscInt *cols;
IS isrows,iscols;
PetscBool flg = PETSC_TRUE;
PetscScalar* v;
PetscRandom rand;
PetscScalar rval;
int procid,size;
/*-----------------------------------------------------------------------------
* Init petsc
*-----------------------------------------------------------------------------*/
PetscInitialize(&argc,&argv,NULL,NULL);
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
MPI_Comm_size(MPI_COMM_WORLD,&size);
ierr = PetscOptionsGetInt(NULL,"-tb",&tb, NULL); CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-tt",&tt, NULL); CHKERRQ(ierr);
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
/*-----------------------------------------------------------------------------
* Set up Matrix MM
*-----------------------------------------------------------------------------*/
ierr = MatCreate(PETSC_COMM_WORLD,&MM);CHKERRQ(ierr);
ierr = MatSetSizes(MM,PETSC_DECIDE,PETSC_DECIDE,tt,tb);CHKERRQ(ierr);
ierr = MatSetType(MM,MATELEMENTAL);CHKERRQ(ierr);
ierr = MatSetFromOptions(MM);CHKERRQ(ierr);
ierr = MatSetUp(MM);CHKERRQ(ierr);
/*-----------------------------------------------------------------------------
* Set Local Values to Matrix MM
*-----------------------------------------------------------------------------*/
ierr=MatGetOwnershipIS(MM,&isrows,&iscols);
ierr=ISGetLocalSize(isrows,&nrows);
ierr=ISGetIndices(isrows,&rows);
ierr=ISGetLocalSize(iscols,&ncols);
ierr=ISGetIndices(iscols,&cols);
ierr=PetscMalloc1(nrows*ncols,&v);
for (i=0; i<nrows; i++) {
for (j=0; j<ncols; j++) {
PetscRandomGetValue(rand,&rval);
v[i*ncols+j] = rval; // PetscReal((i==j)?1:0);
}
}
ierr = MatSetValues(MM,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatAssemblyBegin(MM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(MM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr);
ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr);
ierr = ISDestroy(&isrows);CHKERRQ(ierr);
ierr = ISDestroy(&iscols);CHKERRQ(ierr);
if (flg) {
// ierr=MatView(MM,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/*-----------------------------------------------------------------------------
* Set up sol and rhs vectors s.t rhs = MM*ex_sol
*-----------------------------------------------------------------------------*/
ierr = VecCreate(PETSC_COMM_WORLD,&sol);CHKERRQ(ierr);
ierr = VecSetSizes(sol,PETSC_DECIDE,tb);CHKERRQ(ierr);
ierr = VecSetFromOptions(sol);CHKERRQ(ierr);
ierr = VecDuplicate(sol,&rhs);CHKERRQ(ierr);
ierr = VecDuplicate(sol,&ex_sol);CHKERRQ(ierr);
ierr = VecSet(sol,0.0);CHKERRQ(ierr);
ierr = VecSet(rhs,0.0);CHKERRQ(ierr);
ierr = VecSet(ex_sol,1.0);CHKERRQ(ierr);
ierr = MatMult(MM,ex_sol,rhs);CHKERRQ(ierr);
if(flg){
ierr = VecView(rhs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(sol,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(ex_sol,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/*-----------------------------------------------------------------------------
* Solve system MM*sol = rhs
*-----------------------------------------------------------------------------*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);// optional
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,MM,MM);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-5,1.e-16,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERELEMENTAL);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr);
/*-----------------------------------------------------------------------------
* Check Error
*-----------------------------------------------------------------------------*/
if(flg){
if(procid==0)
std::cout<<"sol is "<<"\n";
MPI_Barrier(MPI_COMM_WORLD);
ierr = VecView(sol,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
ierr = VecAXPY(sol,-1.0,ex_sol);CHKERRQ(ierr);
ierr = VecNorm(sol,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
/*-----------------------------------------------------------------------------
* Destroy objects
*-----------------------------------------------------------------------------*/
ierr = PetscFree(v);CHKERRQ(ierr);
ierr = MatDestroy(&MM);CHKERRQ(ierr);
ierr = VecDestroy(&rhs);CHKERRQ(ierr);
ierr = VecDestroy(&sol);CHKERRQ(ierr);
ierr = VecDestroy(&ex_sol);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
MPI_Barrier(MPI_COMM_WORLD);
PetscFinalize();
return 0;
}
#endif
