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 monitors of type XXXMonitorConverged */
51: struct _n_SlepcConvMon {
52: void *ctx;
53: PetscInt oldnconv; /* previous value of nconv */
54: };
56: /* context for structured eigenproblem matrices created via MatCreateXXX */
57: struct _n_SlepcMatStruct {
58: PetscInt cookie; /* identify which structured matrix */
59: PetscScalar s; /* in BSE sign of the bottom part of the vector */
60: };
61: typedef struct _n_SlepcMatStruct* SlepcMatStruct;
63: #define SLEPC_MAT_STRUCT_BSE 88101
64: #define SLEPC_MAT_STRUCT_HAMILT 88102
66: /*
67: SlepcCheckMatStruct - Check that a given Mat is a structured matrix of the wanted type.
69: Returns true/false in flg if it is given, otherwise yields an error if the check fails.
70: If cookie==0 it will check for any type.
71: */
72: static inline PetscErrorCode SlepcCheckMatStruct(Mat A,PetscInt cookie,PetscBool *flg)
73: {
74: PetscContainer container;
75: SlepcMatStruct mctx;
77: PetscFunctionBegin;
78: if (flg) *flg = PETSC_FALSE;
79: PetscCall(PetscObjectQuery((PetscObject)A,"SlepcMatStruct",(PetscObject*)&container));
80: if (flg && !container) PetscFunctionReturn(PETSC_SUCCESS);
81: PetscCheck(container,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The Mat is not a structured matrix");
82: if (cookie) {
83: PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
84: if (flg && (!mctx || mctx->cookie!=cookie)) PetscFunctionReturn(PETSC_SUCCESS);
85: PetscCheck(mctx && mctx->cookie==cookie,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The type of structured matrix is different from the expected one");
86: }
87: if (flg) *flg = PETSC_TRUE;
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: /*
92: SlepcPrintEigenvalueASCII - Print an eigenvalue on an ASCII viewer.
93: */
94: static inline PetscErrorCode SlepcPrintEigenvalueASCII(PetscViewer viewer,PetscScalar eigr,PetscScalar eigi)
95: {
96: PetscReal re,im;
98: PetscFunctionBegin;
99: #if defined(PETSC_USE_COMPLEX)
100: re = PetscRealPart(eigr);
101: im = PetscImaginaryPart(eigr);
102: (void)eigi;
103: #else
104: re = eigr;
105: im = eigi;
106: #endif
107: /* print zero instead of tiny value */
108: if (PetscAbs(im) && PetscAbs(re)/PetscAbs(im)<PETSC_SMALL) re = 0.0;
109: if (PetscAbs(re) && PetscAbs(im)/PetscAbs(re)<PETSC_SMALL) im = 0.0;
110: /* print as real if imaginary part is zero */
111: if (im!=(PetscReal)0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f%+.5fi",(double)re,(double)im));
112: else PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f",(double)re));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*
117: SlepcViewEigenvector - Outputs an eigenvector xr,xi to a viewer.
118: In complex scalars only xr is written.
119: The name of xr,xi is set before writing, based on the label, the index, and the name of obj.
120: */
121: static inline PetscErrorCode SlepcViewEigenvector(PetscViewer viewer,Vec xr,Vec xi,const char *label,PetscInt index,PetscObject obj)
122: {
123: size_t count;
124: char vname[30];
125: const char *pname;
127: PetscFunctionBegin;
128: PetscCall(PetscObjectGetName(obj,&pname));
129: PetscCall(PetscSNPrintfCount(vname,sizeof(vname),"%s%s",&count,label,PetscDefined(USE_COMPLEX)?"":"r"));
130: count--;
131: PetscCall(PetscSNPrintf(vname+count,sizeof(vname)-count,"%" PetscInt_FMT "_%s",index,pname));
132: PetscCall(PetscObjectSetName((PetscObject)xr,vname));
133: PetscCall(VecView(xr,viewer));
134: #if !defined(PETSC_USE_COMPLEX)
135: vname[count-1] = 'i';
136: PetscCall(PetscObjectSetName((PetscObject)xi,vname));
137: PetscCall(VecView(xi,viewer));
138: #else
139: (void)xi;
140: #endif
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /* Macros for strings with different value in real and complex */
145: #if defined(PETSC_USE_COMPLEX)
146: #define SLEPC_STRING_HERMITIAN "hermitian"
147: #else
148: #define SLEPC_STRING_HERMITIAN "symmetric"
149: #endif
151: /* Private functions that are shared by several classes */
152: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisReference_Private(PetscInt,Vec*,PetscInt*,Vec**);
153: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisDestroy_Private(PetscInt*,Vec**);
154: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcMonitorMakeKey_Internal(const char[],PetscViewerType,PetscViewerFormat,char[]);
155: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscViewerAndFormatCreate_Internal(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**);
157: SLEPC_INTERN PetscErrorCode SlepcCitationsInitialize(void);
158: SLEPC_INTERN PetscErrorCode SlepcInitialize_DynamicLibraries(void);
159: SLEPC_INTERN PetscErrorCode SlepcInitialize_Packages(void);
161: /* Macro to check a sequential Mat (including GPU) */
162: #if !defined(PETSC_USE_DEBUG)
163: #define SlepcMatCheckSeq(h) do {(void)(h);} while (0)
164: #else
165: #if defined(PETSC_HAVE_CUDA)
166: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSECUDA); } while (0)
167: #elif defined(PETSC_HAVE_HIP)
168: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSEHIP); } while (0)
169: #else
170: #define SlepcMatCheckSeq(h) do { PetscCheckTypeName((h),MATSEQDENSE); } while (0)
171: #endif
172: #endif
174: /* Definitions needed to work with GPU kernels */
175: #if defined(PETSC_HAVE_CUPM)
176: #include <petscdevice_cupm.h>
178: #define X_AXIS 0
179: #define Y_AXIS 1
181: #define SLEPC_TILE_SIZE_X 32
182: #define SLEPC_BLOCK_SIZE_X 128
183: #define SLEPC_TILE_SIZE_Y 32
184: #define SLEPC_BLOCK_SIZE_Y 128
186: static inline PetscErrorCode SlepcKernelSetGrid1D(PetscInt rows,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount)
187: {
188: int card;
189: #if defined(PETSC_HAVE_CUDA)
190: struct cudaDeviceProp devprop;
191: #elif defined(PETSC_HAVE_HIP)
192: hipDeviceProp_t devprop;
193: #endif
195: PetscFunctionBegin;
196: #if defined(PETSC_HAVE_CUDA)
197: PetscCallCUDA(cudaGetDevice(&card));
198: PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
199: #elif defined(PETSC_HAVE_HIP)
200: PetscCallHIP(hipGetDevice(&card));
201: PetscCallHIP(hipGetDeviceProperties(&devprop,card));
202: #endif
203: *dimGrid_xcount = 1;
205: /* X axis */
206: dimGrid->x = 1;
207: dimBlock->x = SLEPC_BLOCK_SIZE_X;
208: if (rows>SLEPC_BLOCK_SIZE_X) dimGrid->x = (rows+SLEPC_BLOCK_SIZE_X-1)/SLEPC_BLOCK_SIZE_X;
209: else dimBlock->x = rows;
210: if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
211: *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
212: dimGrid->x = devprop.maxGridSize[X_AXIS];
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static inline PetscErrorCode SlepcKernelSetGrid2DTiles(PetscInt rows,PetscInt cols,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount,PetscInt *dimGrid_ycount)
218: {
219: int card;
220: #if defined(PETSC_HAVE_CUDA)
221: struct cudaDeviceProp devprop;
222: #elif defined(PETSC_HAVE_HIP)
223: hipDeviceProp_t devprop;
224: #endif
226: PetscFunctionBegin;
227: #if defined(PETSC_HAVE_CUDA)
228: PetscCallCUDA(cudaGetDevice(&card));
229: PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
230: #elif defined(PETSC_HAVE_HIP)
231: PetscCallHIP(hipGetDevice(&card));
232: PetscCallHIP(hipGetDeviceProperties(&devprop,card));
233: #endif
234: *dimGrid_xcount = *dimGrid_ycount = 1;
236: /* X axis */
237: dimGrid->x = 1;
238: dimBlock->x = SLEPC_BLOCK_SIZE_X;
239: 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);
240: else dimBlock->x = (rows+SLEPC_TILE_SIZE_X-1)/SLEPC_TILE_SIZE_X;
241: if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
242: *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
243: dimGrid->x = devprop.maxGridSize[X_AXIS];
244: }
246: /* Y axis */
247: dimGrid->y = 1;
248: dimBlock->y = SLEPC_BLOCK_SIZE_Y;
249: 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);
250: else dimBlock->y = (cols+SLEPC_TILE_SIZE_Y-1)/SLEPC_TILE_SIZE_Y;
251: if (dimGrid->y>(unsigned)devprop.maxGridSize[Y_AXIS]) {
252: *dimGrid_ycount = (dimGrid->y+(devprop.maxGridSize[Y_AXIS]-1))/devprop.maxGridSize[Y_AXIS];
253: dimGrid->y = devprop.maxGridSize[Y_AXIS];
254: }
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
257: #undef X_AXIS
258: #undef Y_AXIS
259: #endif