Actual source code: slepcimpl.h
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: #pragma once
13: #include <slepcsys.h>
14: #include <petsc/private/petscimpl.h>
16: /* SUBMANSEC = Sys */
18: SLEPC_INTERN PetscBool SlepcBeganPetsc;
20: /* SlepcSwap - swap two variables a,b of the same type using a temporary variable t */
21: #define SlepcSwap(a,b,t) do {t=a;a=b;b=t;} while (0)
23: /*MC
24: SlepcHeaderCreate - Creates a SLEPc object
26: Input Parameters:
27: + classid - the classid associated with this object
28: . class_name - string name of class; should be static
29: . descr - string containing short description; should be static
30: . mansec - string indicating section in manual pages; should be static
31: . comm - the MPI Communicator
32: . destroy - the destroy routine for this object
33: - view - the view routine for this object
35: Output Parameter:
36: . h - the newly created object
38: Note:
39: This is equivalent to `PetscHeaderCreate()` but makes sure that `SlepcInitialize()`
40: has been called.
42: Level: developer
43: M*/
44: #define SlepcHeaderCreate(h,classid,class_name,descr,mansec,comm,destroy,view) \
45: ((PetscErrorCode)((!SlepcInitializeCalled && \
46: PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,PETSC_ERR_ORDER,PETSC_ERROR_INITIAL, \
47: "Must call SlepcInitialize instead of PetscInitialize to use SLEPc classes")) || \
48: PetscHeaderCreate(h,classid,class_name,descr,mansec,comm,destroy,view)))
50: /* context for structured eigenproblem matrices created via MatCreateXXX */
51: struct _n_SlepcMatStruct {
52: PetscInt cookie; /* identify which structured matrix */
53: PetscScalar s; /* in BSE sign of the bottom part of the vector */
54: };
55: typedef struct _n_SlepcMatStruct* SlepcMatStruct;
57: #define SLEPC_MAT_STRUCT_BSE 88101
58: #define SLEPC_MAT_STRUCT_HAMILT 88102
59: #define SLEPC_MAT_STRUCT_LREP 88103
61: /*
62: SlepcCheckMatStruct - Check that a given Mat is a structured matrix of the wanted type.
64: Returns true/false in flg if it is given, otherwise yields an error if the check fails.
65: If cookie==0 it will check for any type.
66: */
67: static inline PetscErrorCode SlepcCheckMatStruct(Mat A,PetscInt cookie,PetscBool *flg)
68: {
69: PetscContainer container;
70: SlepcMatStruct mctx;
72: PetscFunctionBegin;
73: if (flg) *flg = PETSC_FALSE;
74: PetscCall(PetscObjectQuery((PetscObject)A,"SlepcMatStruct",(PetscObject*)&container));
75: if (flg && !container) PetscFunctionReturn(PETSC_SUCCESS);
76: PetscCheck(container,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The Mat is not a structured matrix");
77: if (cookie) {
78: PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
79: if (flg && (!mctx || mctx->cookie!=cookie)) PetscFunctionReturn(PETSC_SUCCESS);
80: PetscCheck(mctx && mctx->cookie==cookie,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The type of structured matrix is different from the expected one");
81: }
82: if (flg) *flg = PETSC_TRUE;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*
87: SlepcCheckMatLREPReduced - Check whether a given LREP Mat is of reduced form or not.
89: True if the matrix has reduced form [0 K; M 0], otherwise it has the form [A B; -B -A].
90: */
91: static inline PetscErrorCode SlepcCheckMatLREPReduced(Mat A,PetscBool *flg)
92: {
93: Mat A00;
95: PetscFunctionBegin;
96: PetscAssertPointer(flg,2);
97: *flg = PETSC_FALSE;
98: PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_LREP,NULL)); /* error out if not LREP */
99: PetscCall(MatNestGetSubMat(A,0,0,&A00));
100: if (A00==NULL) *flg = PETSC_TRUE;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*
105: SlepcPrintEigenvalueASCII - Print an eigenvalue on an ASCII viewer.
106: */
107: static inline PetscErrorCode SlepcPrintEigenvalueASCII(PetscViewer viewer,PetscScalar eigr,PetscScalar eigi)
108: {
109: PetscReal re,im;
111: PetscFunctionBegin;
112: #if defined(PETSC_USE_COMPLEX)
113: re = PetscRealPart(eigr);
114: im = PetscImaginaryPart(eigr);
115: (void)eigi;
116: #else
117: re = eigr;
118: im = eigi;
119: #endif
120: /* print zero instead of tiny value */
121: if (PetscAbs(im) && PetscAbs(re)/PetscAbs(im)<PETSC_SMALL) re = 0.0;
122: if (PetscAbs(re) && PetscAbs(im)/PetscAbs(re)<PETSC_SMALL) im = 0.0;
123: /* print as real if imaginary part is zero */
124: if (im!=(PetscReal)0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f%+.5fi",(double)re,(double)im));
125: else PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f",(double)re));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*
130: SlepcViewEigenvector - Outputs an eigenvector xr,xi to a viewer.
131: In complex scalars only xr is written.
132: The name of xr,xi is set before writing, based on the label, the index, and the name of obj.
133: */
134: static inline PetscErrorCode SlepcViewEigenvector(PetscViewer viewer,Vec xr,Vec xi,const char *label,PetscInt index,PetscObject obj)
135: {
136: size_t count;
137: char vname[30];
138: const char *pname;
140: PetscFunctionBegin;
141: PetscCall(PetscObjectGetName(obj,&pname));
142: PetscCall(PetscSNPrintfCount(vname,sizeof(vname),"%s%s",&count,label,PetscDefined(USE_COMPLEX)?"":"r"));
143: count--;
144: PetscCall(PetscSNPrintf(vname+count,sizeof(vname)-count,"%" PetscInt_FMT "_%s",index,pname));
145: PetscCall(PetscObjectSetName((PetscObject)xr,vname));
146: PetscCall(VecView(xr,viewer));
147: if (!PetscDefined(USE_COMPLEX)) {
148: vname[count-1] = 'i';
149: PetscCall(PetscObjectSetName((PetscObject)xi,vname));
150: PetscCall(VecView(xi,viewer));
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /* Macros for strings with different value in real and complex */
156: #if defined(PETSC_USE_COMPLEX)
157: #define SLEPC_STRING_HERMITIAN "hermitian"
158: #else
159: #define SLEPC_STRING_HERMITIAN "symmetric"
160: #endif
162: /* Private functions that are shared by several classes */
163: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisReference_Private(PetscInt,Vec*,PetscInt*,Vec**);
164: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisDestroy_Private(PetscInt*,Vec**);
165: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcMonitorMakeKey_Internal(const char[],PetscViewerType,PetscViewerFormat,char[]);
166: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscViewerAndFormatCreate_Internal(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**);
168: SLEPC_INTERN PetscErrorCode SlepcCitationsInitialize(void);
169: SLEPC_INTERN PetscErrorCode SlepcInitialize_DynamicLibraries(void);
170: SLEPC_INTERN PetscErrorCode SlepcInitialize_Packages(void);
172: /* Macro to check a sequential Mat (including GPU) */
173: #if !defined(PETSC_USE_DEBUG)
174: #define SlepcMatCheckSeq(h) do {(void)(h);} while (0)
175: #else
176: #if defined(PETSC_HAVE_CUDA)
177: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSECUDA); } while (0)
178: #elif defined(PETSC_HAVE_HIP)
179: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSEHIP); } while (0)
180: #else
181: #define SlepcMatCheckSeq(h) do { PetscCheckTypeName((h),MATSEQDENSE); } while (0)
182: #endif
183: #endif
185: /* Definitions needed to work with GPU kernels */
186: #if defined(PETSC_HAVE_CUPM)
187: #include <petscdevice_cupm.h>
189: #define X_AXIS 0
190: #define Y_AXIS 1
192: #define SLEPC_TILE_SIZE_X 32
193: #define SLEPC_BLOCK_SIZE_X 128
194: #define SLEPC_TILE_SIZE_Y 32
195: #define SLEPC_BLOCK_SIZE_Y 128
197: static inline PetscErrorCode SlepcKernelSetGrid1D(PetscInt rows,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount)
198: {
199: int card;
200: #if defined(PETSC_HAVE_CUDA)
201: struct cudaDeviceProp devprop;
202: #elif defined(PETSC_HAVE_HIP)
203: hipDeviceProp_t devprop;
204: #endif
206: PetscFunctionBegin;
207: #if defined(PETSC_HAVE_CUDA)
208: PetscCallCUDA(cudaGetDevice(&card));
209: PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
210: #elif defined(PETSC_HAVE_HIP)
211: PetscCallHIP(hipGetDevice(&card));
212: PetscCallHIP(hipGetDeviceProperties(&devprop,card));
213: #endif
214: *dimGrid_xcount = 1;
216: /* X axis */
217: dimGrid->x = 1;
218: dimBlock->x = SLEPC_BLOCK_SIZE_X;
219: if (rows>SLEPC_BLOCK_SIZE_X) dimGrid->x = (rows+SLEPC_BLOCK_SIZE_X-1)/SLEPC_BLOCK_SIZE_X;
220: else dimBlock->x = rows;
221: if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
222: *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
223: dimGrid->x = devprop.maxGridSize[X_AXIS];
224: }
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static inline PetscErrorCode SlepcKernelSetGrid2DTiles(PetscInt rows,PetscInt cols,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount,PetscInt *dimGrid_ycount)
229: {
230: int card;
231: #if defined(PETSC_HAVE_CUDA)
232: struct cudaDeviceProp devprop;
233: #elif defined(PETSC_HAVE_HIP)
234: hipDeviceProp_t devprop;
235: #endif
237: PetscFunctionBegin;
238: #if defined(PETSC_HAVE_CUDA)
239: PetscCallCUDA(cudaGetDevice(&card));
240: PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
241: #elif defined(PETSC_HAVE_HIP)
242: PetscCallHIP(hipGetDevice(&card));
243: PetscCallHIP(hipGetDeviceProperties(&devprop,card));
244: #endif
245: *dimGrid_xcount = *dimGrid_ycount = 1;
247: /* X axis */
248: dimGrid->x = 1;
249: dimBlock->x = SLEPC_BLOCK_SIZE_X;
250: if (rows>SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X) dimGrid->x = (rows+SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X-1)/(SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X);
251: else dimBlock->x = (rows+SLEPC_TILE_SIZE_X-1)/SLEPC_TILE_SIZE_X;
252: if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
253: *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
254: dimGrid->x = devprop.maxGridSize[X_AXIS];
255: }
257: /* Y axis */
258: dimGrid->y = 1;
259: dimBlock->y = SLEPC_BLOCK_SIZE_Y;
260: if (cols>SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y) dimGrid->y = (cols+SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y-1)/(SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y);
261: else dimBlock->y = (cols+SLEPC_TILE_SIZE_Y-1)/SLEPC_TILE_SIZE_Y;
262: if (dimGrid->y>(unsigned)devprop.maxGridSize[Y_AXIS]) {
263: *dimGrid_ycount = (dimGrid->y+(devprop.maxGridSize[Y_AXIS]-1))/devprop.maxGridSize[Y_AXIS];
264: dimGrid->y = devprop.maxGridSize[Y_AXIS];
265: }
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
268: #undef X_AXIS
269: #undef Y_AXIS
270: #endif