Actual source code: test2.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Tests MatNormEstimate() on the Brusselator matrix.\n\n"
     12:   "The command line options are:\n"
     13:   "  -n <n>, where <n> = block dimension of the 2x2 block matrix.\n"
     14:   "  -L <L>, where <L> = bifurcation parameter.\n"
     15:   "  -alpha <alpha>, -beta <beta>, -delta1 <delta1>,  -delta2 <delta2>,\n"
     16:   "       where <alpha> <beta> <delta1> <delta2> = model parameters.\n\n";

     18: #include <slepcsys.h>

     20: /*
     21:    The Brusselator matrix is

     23:         A = [ tau1*T+(beta-1)*I     alpha^2*I
     24:                   -beta*I        tau2*T-alpha^2*I ],

     26:    where

     28:         T = tridiag{1,-2,1}
     29:         h = 1/(n+1)
     30:         tau1 = delta1/(h*L)^2
     31:         tau2 = delta2/(h*L)^2
     32:  */

     34: int main(int argc,char **argv)
     35: {
     36:   Mat            A,T1,T2,D1,D2,mats[4];
     37:   PetscScalar    alpha,beta,tau1,tau2,delta1,delta2,L,h;
     38:   PetscReal      nrm;
     39:   PetscInt       N=30,i,Istart,Iend;

     41:   PetscFunctionBeginUser;
     42:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     43:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));

     45:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     46:         Generate the matrix
     47:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     49:   alpha  = 2.0;
     50:   beta   = 5.45;
     51:   delta1 = 0.008;
     52:   delta2 = 0.004;
     53:   L      = 0.51302;

     55:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
     56:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
     57:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL));
     58:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
     59:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));

     61:   h    = 1.0 / (PetscReal)(N+1);
     62:   tau1 = delta1 / ((h*L)*(h*L));
     63:   tau2 = delta2 / ((h*L)*(h*L));

     65:   /* Create matrices T1, T2 */
     66:   PetscCall(MatCreate(PETSC_COMM_WORLD,&T1));
     67:   PetscCall(MatSetSizes(T1,PETSC_DECIDE,PETSC_DECIDE,N,N));
     68:   PetscCall(MatSetFromOptions(T1));
     69:   PetscCall(MatGetOwnershipRange(T1,&Istart,&Iend));
     70:   for (i=Istart;i<Iend;i++) {
     71:     if (i>0) PetscCall(MatSetValue(T1,i,i-1,1.0,INSERT_VALUES));
     72:     if (i<N-1) PetscCall(MatSetValue(T1,i,i+1,1.0,INSERT_VALUES));
     73:     PetscCall(MatSetValue(T1,i,i,-2.0,INSERT_VALUES));
     74:   }
     75:   PetscCall(MatAssemblyBegin(T1,MAT_FINAL_ASSEMBLY));
     76:   PetscCall(MatAssemblyEnd(T1,MAT_FINAL_ASSEMBLY));

     78:   PetscCall(MatDuplicate(T1,MAT_COPY_VALUES,&T2));
     79:   PetscCall(MatScale(T1,tau1));
     80:   PetscCall(MatShift(T1,beta-1.0));
     81:   PetscCall(MatScale(T2,tau2));
     82:   PetscCall(MatShift(T2,-alpha*alpha));

     84:   /* Create matrices D1, D2 */
     85:   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,alpha*alpha,&D1));
     86:   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,-beta,&D2));

     88:   /* Create the nest matrix */
     89:   mats[0] = T1;
     90:   mats[1] = D1;
     91:   mats[2] = D2;
     92:   mats[3] = T2;
     93:   PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,mats,&A));

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:         Estimate the norm
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   PetscCall(MatNormEstimate(A,NULL,NULL,&nrm));
    100:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator matrix, n=%" PetscInt_FMT " - estimated norm = %g\n\n",N,(double)nrm));

    102:   PetscCall(MatDestroy(&A));
    103:   PetscCall(MatDestroy(&T1));
    104:   PetscCall(MatDestroy(&T2));
    105:   PetscCall(MatDestroy(&D1));
    106:   PetscCall(MatDestroy(&D2));
    107:   PetscCall(SlepcFinalize());
    108:   return 0;
    109: }

    111: /*TEST

    113:    test:
    114:      suffix: 1
    115:      requires: !complex

    117:    test:
    118:      suffix: 1_complex
    119:      requires: complex

    121: TEST*/