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*/