Actual source code: test18.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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[] = "Test DSSynchronize() on a NHEP.\n\n";
13: #include <slepcds.h>
15: PetscErrorCode CheckArray(PetscScalar *A,const char *label,PetscInt k)
16: {
18: PetscInt j;
19: PetscMPIInt p,size,rank;
20: PetscScalar dif,*buf;
21: PetscReal error;
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: if (rank) {
27: MPI_Send(A,k,MPIU_SCALAR,0,111,PETSC_COMM_WORLD);
28: } else {
29: PetscMalloc1(k,&buf);
30: for (p=1;p<size;p++) {
31: MPI_Recv(buf,k,MPIU_SCALAR,p,111,PETSC_COMM_WORLD,MPI_STATUS_IGNORE);
32: dif = 0.0;
33: for (j=0;j<k;j++) dif += A[j]-buf[j];
34: error = PetscAbsScalar(dif);
35: if (error>10*PETSC_MACHINE_EPSILON) {
36: PetscPrintf(PETSC_COMM_WORLD," Array %s differs in proc %d: %g\n",label,(int)p,(double)error);
37: }
38: }
39: PetscFree(buf);
40: }
41: return(0);
42: }
44: int main(int argc,char **argv)
45: {
47: DS ds;
48: SlepcSC sc;
49: PetscScalar *A,*Q,*wr,*wi;
50: PetscReal re,im;
51: PetscInt i,j,n=10;
52: PetscMPIInt size;
54: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
55: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
56: PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type NHEP - dimension %D.\n",n);
58: /* Create DS object */
59: DSCreate(PETSC_COMM_WORLD,&ds);
60: DSSetType(ds,DSNHEP);
61: DSSetFromOptions(ds);
62: DSAllocate(ds,n);
63: DSSetDimensions(ds,n,0,0,0);
65: /* Fill with Grcar matrix */
66: DSGetArray(ds,DS_MAT_A,&A);
67: for (i=1;i<n;i++) A[i+(i-1)*n]=-1.0;
68: for (j=0;j<4;j++) {
69: for (i=0;i<n-j;i++) A[i+(i+j)*n]=1.0;
70: }
71: DSRestoreArray(ds,DS_MAT_A,&A);
72: DSSetState(ds,DS_STATE_INTERMEDIATE);
74: /* Solve */
75: PetscMalloc2(n,&wr,n,&wi);
76: DSGetSlepcSC(ds,&sc);
77: sc->comparison = SlepcCompareLargestMagnitude;
78: sc->comparisonctx = NULL;
79: sc->map = NULL;
80: sc->mapobj = NULL;
81: DSSolve(ds,wr,wi);
82: DSSort(ds,wr,wi,NULL,NULL,NULL);
84: /* Print eigenvalues */
85: PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues =\n");
86: for (i=0;i<n;i++) {
87: #if defined(PETSC_USE_COMPLEX)
88: re = PetscRealPart(wr[i]);
89: im = PetscImaginaryPart(wr[i]);
90: #else
91: re = wr[i];
92: im = wi[i];
93: #endif
94: if (PetscAbs(im)<1e-10) {
95: PetscPrintf(PETSC_COMM_WORLD," %.5f\n",(double)re);
96: } else {
97: PetscPrintf(PETSC_COMM_WORLD," %.5f%+.5fi\n",(double)re,(double)im);
98: }
99: }
101: /* Synchronize data and check */
102: DSSynchronize(ds,wr,wi);
103: MPI_Comm_size(PETSC_COMM_WORLD,&size);
104: if (size>1) {
105: CheckArray(wr,"wr",n);
106: #if !defined(PETSC_USE_COMPLEX)
107: CheckArray(wi,"wi",n);
108: #endif
109: DSGetArray(ds,DS_MAT_A,&A);
110: CheckArray(A,"A",n*n);
111: DSRestoreArray(ds,DS_MAT_A,&A);
112: DSGetArray(ds,DS_MAT_Q,&Q);
113: CheckArray(Q,"Q",n*n);
114: DSRestoreArray(ds,DS_MAT_Q,&Q);
115: }
117: PetscFree2(wr,wi);
118: DSDestroy(&ds);
119: SlepcFinalize();
120: return ierr;
121: }
123: /*TEST
125: test:
126: suffix: 1
127: nsize: {{1 2 3}}
128: args: -ds_parallel redundant
129: requires: !complex
130: output_file: output/test18_1.out
132: test:
133: suffix: 2
134: nsize: {{1 2 3}}
135: args: -ds_parallel synchronized
136: requires: !complex
137: output_file: output/test18_1.out
139: TEST*/