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