static char help[] = "Solves a tridiagonal linear system with KSP.\n\n"; /*T Concepts: KSP^solving a system of linear equations Processors: 1 T*/ /* Include "petscksp.h" so that we can use KSP solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners Note: The corresponding parallel example is ex23.c */ #include #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc,char **args) { Vec x, b, u; /* approx solution, RHS, exact solution */ Mat A; /* linear system matrix */ KSP ksp; /* linear solver context */ PC pc; /* preconditioner context */ PetscReal norm,tol=1.e-14; /* norm of solution error */ PetscErrorCode ierr; PetscInt i,n = 10,col[3],its; PetscMPIInt size; PetscScalar neg_one = -1.0,one = 1.0,value[3]; PetscBool nonzeroguess = PETSC_FALSE; // profiling structs PetscClassId user_classid; PetscLogStage stage_build, stage_solve; PetscLogEvent event_build, event_solve; PetscEventPerfInfo info_build, info_solve; PetscStageLog stageLog; PetscInitialize(&argc,&args,(char*)0,help); ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!"); ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetBool(NULL,"-nonzero_guess",&nonzeroguess,NULL);CHKERRQ(ierr); // registrations PetscClassIdRegister("ex1", &user_classid); PetscLogStageRegister("Assembly", &stage_build); PetscLogStageRegister("Solution", &stage_solve); PetscLogEventRegister("build", user_classid, &event_build); PetscLogEventRegister("solve", user_classid, &event_solve); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute the matrix and right-hand-side vector that define the linear system, Ax = b. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create vectors. Note that we form 1 vector from scratch and then duplicate as needed. */ ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr); ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr); ierr = VecSetFromOptions(x);CHKERRQ(ierr); ierr = VecDuplicate(x,&b);CHKERRQ(ierr); ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* Create matrix. When using MatCreate(), the matrix format can be specified at runtime. Performance tuning note: For problems of substantial size, preallocation of matrix memory is crucial for attaining good performance. See the matrix chapter of the users manual for details. */ 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); /* Assemble matrix */ PetscLogStagePush(stage_build); PetscLogEventBegin(event_build, 0, 0, 0, 0); value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; for (i=1; i -pc_type -ksp_monitor -ksp_rtol These options will override those specified above as long as KSPSetFromOptions() is called _after_ any other customization routines. */ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); if (nonzeroguess) { PetscScalar p = .5; ierr = VecSet(x,p);CHKERRQ(ierr); ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Solve linear system */ ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); /* View solver info; we could instead use the option -ksp_view to print this info to the screen at the conclusion of KSPSolve(). */ ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Check solution and clean up - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Check the error */ ierr = VecAXPY(x,neg_one,u);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); if (norm > tol) { ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); } /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */ ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = KSPDestroy(&ksp);CHKERRQ(ierr); PetscLogEventEnd(event_solve, 0, 0, 0, 0); PetscLogStagePop(); PetscPrintf(PETSC_COMM_WORLD, "profiling information\n"); PetscPrintf(PETSC_COMM_WORLD, "=====================\n\n"); // these infos are empty PetscLogEventGetPerfInfo(stage_build, event_build, &info_build); PetscPrintf(PETSC_COMM_WORLD, "info_build.time = %g\n", info_build.time); PetscLogEventGetPerfInfo(stage_solve, event_solve, &info_solve); PetscPrintf(PETSC_COMM_WORLD, "info_solve.time = %g\n", info_solve.time); // however these instructions work PetscLogGetStageLog(&stageLog); PetscPrintf(PETSC_COMM_WORLD, "stageLog->stageInfo[stage_build].perfInfo.time = %g\n", stageLog->stageInfo[stage_build].perfInfo.time); PetscPrintf(PETSC_COMM_WORLD, "stageLog->stageInfo[stage_solve].perfInfo.time = %g\n", stageLog->stageInfo[stage_solve].perfInfo.time); /* Always call PetscFinalize() before exiting a program. This routine - finalizes the PETSc libraries as well as MPI - provides summary and diagnostic information if certain runtime options are chosen (e.g., -log_summary). */ ierr = PetscFinalize(); return 0; }