Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test DSSynchronize() on a NHEP.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 40 : PetscErrorCode CheckArray(PetscScalar *A,const char *label,PetscInt k)
16 : {
17 40 : PetscInt j;
18 40 : PetscMPIInt p,size,rank;
19 40 : PetscScalar dif,*buf;
20 40 : PetscReal error;
21 :
22 40 : PetscFunctionBeginUser;
23 40 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
24 40 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
25 40 : if (rank) PetscCallMPI(MPI_Send(A,k,MPIU_SCALAR,0,111,PETSC_COMM_WORLD));
26 : else {
27 16 : PetscCall(PetscMalloc1(k,&buf));
28 40 : for (p=1;p<size;p++) {
29 24 : PetscCallMPI(MPI_Recv(buf,k,MPIU_SCALAR,p,111,PETSC_COMM_WORLD,MPI_STATUS_IGNORE));
30 : dif = 0.0;
31 1344 : for (j=0;j<k;j++) dif += A[j]-buf[j];
32 24 : error = PetscAbsScalar(dif);
33 24 : if (error>10*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Array %s differs in proc %d: %g\n",label,(int)p,(double)error));
34 : }
35 16 : PetscCall(PetscFree(buf));
36 : }
37 40 : PetscFunctionReturn(PETSC_SUCCESS);
38 : }
39 :
40 12 : int main(int argc,char **argv)
41 : {
42 12 : DS ds;
43 12 : SlepcSC sc;
44 12 : PetscScalar *A,*Q,*wr,*wi;
45 12 : PetscReal re,im;
46 12 : PetscInt i,j,n=10;
47 12 : PetscMPIInt size;
48 :
49 12 : PetscFunctionBeginUser;
50 12 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
51 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
52 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type NHEP - dimension %" PetscInt_FMT ".\n",n));
53 :
54 : /* Create DS object */
55 12 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
56 12 : PetscCall(DSSetType(ds,DSNHEP));
57 12 : PetscCall(DSSetFromOptions(ds));
58 12 : PetscCall(DSAllocate(ds,n));
59 12 : PetscCall(DSSetDimensions(ds,n,0,0));
60 :
61 : /* Fill with Grcar matrix */
62 12 : PetscCall(DSGetArray(ds,DS_MAT_A,&A));
63 120 : for (i=1;i<n;i++) A[i+(i-1)*n]=-1.0;
64 60 : for (j=0;j<4;j++) {
65 456 : for (i=0;i<n-j;i++) A[i+(i+j)*n]=1.0;
66 : }
67 12 : PetscCall(DSRestoreArray(ds,DS_MAT_A,&A));
68 12 : PetscCall(DSSetState(ds,DS_STATE_INTERMEDIATE));
69 :
70 : /* Solve */
71 12 : PetscCall(PetscMalloc2(n,&wr,n,&wi));
72 12 : PetscCall(DSGetSlepcSC(ds,&sc));
73 12 : sc->comparison = SlepcCompareLargestMagnitude;
74 12 : sc->comparisonctx = NULL;
75 12 : sc->map = NULL;
76 12 : sc->mapobj = NULL;
77 12 : PetscCall(DSSolve(ds,wr,wi));
78 12 : PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));
79 :
80 : /* Print eigenvalues */
81 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues =\n"));
82 132 : for (i=0;i<n;i++) {
83 : #if defined(PETSC_USE_COMPLEX)
84 : re = PetscRealPart(wr[i]);
85 : im = PetscImaginaryPart(wr[i]);
86 : #else
87 120 : re = wr[i];
88 120 : im = wi[i];
89 : #endif
90 120 : if (PetscAbs(im)<1e-10) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %.5f\n",(double)re));
91 120 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %.5f%+.5fi\n",(double)re,(double)im));
92 : }
93 :
94 : /* Synchronize data and check */
95 12 : PetscCall(DSSynchronize(ds,wr,wi));
96 12 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
97 12 : if (size>1) {
98 10 : PetscCall(CheckArray(wr,"wr",n));
99 : #if !defined(PETSC_USE_COMPLEX)
100 10 : PetscCall(CheckArray(wi,"wi",n));
101 : #endif
102 10 : PetscCall(DSGetArray(ds,DS_MAT_A,&A));
103 10 : PetscCall(CheckArray(A,"A",n*n));
104 10 : PetscCall(DSRestoreArray(ds,DS_MAT_A,&A));
105 10 : PetscCall(DSGetArray(ds,DS_MAT_Q,&Q));
106 10 : PetscCall(CheckArray(Q,"Q",n*n));
107 10 : PetscCall(DSRestoreArray(ds,DS_MAT_Q,&Q));
108 : }
109 :
110 12 : PetscCall(PetscFree2(wr,wi));
111 12 : PetscCall(DSDestroy(&ds));
112 12 : PetscCall(SlepcFinalize());
113 : return 0;
114 : }
115 :
116 : /*TEST
117 :
118 : testset:
119 : nsize: {{1 2 3}}
120 : filter: sed -e "s/[+-]\([0-9]\.[0-9]*i\)/+-\\1/" | sed -e "s/1.58254/1.58255/" | sed -e "s/1.75989/1.75988/" | sed -e "s/1.01945/1.01946/"
121 : output_file: output/test18_1.out
122 : test:
123 : suffix: 1
124 : args: -ds_parallel redundant
125 : test:
126 : suffix: 2
127 : args: -ds_parallel synchronized
128 :
129 : TEST*/
|
