1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17from _ceed_cffi import ffi, lib 18import tempfile 19import numpy as np 20import contextlib 21from .ceed_constants import MEM_HOST, COPY_VALUES, NORM_2 22 23# ------------------------------------------------------------------------------ 24 25 26class Vector(): 27 """Ceed Vector: storing and manipulating vectors.""" 28 29 # Attributes 30 _ceed = ffi.NULL 31 _pointer = ffi.NULL 32 33 # Constructor 34 def __init__(self, ceed, size): 35 # CeedVector object 36 self._pointer = ffi.new("CeedVector *") 37 38 # Reference to Ceed 39 self._ceed = ceed 40 41 # libCEED call 42 err_code = lib.CeedVectorCreate( 43 self._ceed._pointer[0], size, self._pointer) 44 self._ceed._check_error(err_code) 45 46 # Destructor 47 def __del__(self): 48 # libCEED call 49 err_code = lib.CeedVectorDestroy(self._pointer) 50 self._ceed._check_error(err_code) 51 52 # Representation 53 def __repr__(self): 54 return "<CeedVector instance at " + hex(id(self)) + ">" 55 56 # String conversion for print() to stdout 57 def __str__(self): 58 """View a Vector via print().""" 59 60 # libCEED call 61 fmt = ffi.new("char[]", "%f".encode('ascii')) 62 with tempfile.NamedTemporaryFile() as key_file: 63 with open(key_file.name, 'r+') as stream_file: 64 stream = ffi.cast("FILE *", stream_file) 65 66 err_code = lib.CeedVectorView(self._pointer[0], fmt, stream) 67 self._ceed._check_error(err_code) 68 69 stream_file.seek(0) 70 out_string = stream_file.read() 71 72 return out_string 73 74 # Set Vector's data array 75 def set_array(self, array, memtype=MEM_HOST, cmode=COPY_VALUES): 76 """Set the array used by a Vector, freeing any previously allocated 77 array if applicable. 78 79 Args: 80 *array: Numpy or Numba array to be used 81 **memtype: memory type of the array being passed, default CEED_MEM_HOST 82 **cmode: copy mode for the array, default CEED_COPY_VALUES""" 83 84 # Setup the numpy array for the libCEED call 85 array_pointer = ffi.new("CeedScalar *") 86 if memtype == MEM_HOST: 87 array_pointer = ffi.cast( 88 "CeedScalar *", 89 array.__array_interface__['data'][0]) 90 else: 91 array_pointer = ffi.cast( 92 "CeedScalar *", 93 array.__cuda_array_interface__['data'][0]) 94 95 # libCEED call 96 err_code = lib.CeedVectorSetArray( 97 self._pointer[0], memtype, cmode, array_pointer) 98 self._ceed._check_error(err_code) 99 100 # Get Vector's data array 101 def get_array(self, memtype=MEM_HOST): 102 """Get read/write access to a Vector via the specified memory type. 103 104 Args: 105 **memtype: memory type of the array being passed, default CEED_MEM_HOST 106 107 Returns: 108 *array: Numpy or Numba array""" 109 110 # Retrieve the length of the array 111 length_pointer = ffi.new("CeedInt *") 112 err_code = lib.CeedVectorGetLength(self._pointer[0], length_pointer) 113 self._ceed._check_error(err_code) 114 115 # Setup the pointer's pointer 116 array_pointer = ffi.new("CeedScalar **") 117 118 # libCEED call 119 err_code = lib.CeedVectorGetArray( 120 self._pointer[0], memtype, array_pointer) 121 self._ceed._check_error(err_code) 122 123 # Return array created from buffer 124 if memtype == MEM_HOST: 125 # Create buffer object from returned pointer 126 buff = ffi.buffer( 127 array_pointer[0], 128 ffi.sizeof("CeedScalar") * 129 length_pointer[0]) 130 # return Numpy array 131 return np.frombuffer(buff, dtype="float64") 132 else: 133 # CUDA array interface 134 # https://numba.pydata.org/numba-doc/latest/cuda/cuda_array_interface.html 135 import numba.cuda as nbcuda 136 desc = { 137 'shape': (length_pointer[0]), 138 'typestr': '>f8', 139 'data': (int(ffi.cast("intptr_t", array_pointer[0])), False), 140 'version': 2 141 } 142 # return Numba array 143 return nbcuda.from_cuda_array_interface(desc) 144 145 # Get Vector's data array in read-only mode 146 def get_array_read(self, memtype=MEM_HOST): 147 """Get read-only access to a Vector via the specified memory type. 148 149 Args: 150 **memtype: memory type of the array being passed, default CEED_MEM_HOST 151 152 Returns: 153 *array: Numpy or Numba array""" 154 155 # Retrieve the length of the array 156 length_pointer = ffi.new("CeedInt *") 157 err_code = lib.CeedVectorGetLength(self._pointer[0], length_pointer) 158 self._ceed._check_error(err_code) 159 160 # Setup the pointer's pointer 161 array_pointer = ffi.new("CeedScalar **") 162 163 # libCEED call 164 err_code = lib.CeedVectorGetArrayRead( 165 self._pointer[0], memtype, array_pointer) 166 self._ceed._check_error(err_code) 167 168 # Return array created from buffer 169 if memtype == MEM_HOST: 170 # Create buffer object from returned pointer 171 buff = ffi.buffer( 172 array_pointer[0], 173 ffi.sizeof("CeedScalar") * 174 length_pointer[0]) 175 # return read only Numpy array 176 ret = np.frombuffer(buff, dtype="float64") 177 ret.flags['WRITEABLE'] = False 178 return ret 179 else: 180 # CUDA array interface 181 # https://numba.pydata.org/numba-doc/latest/cuda/cuda_array_interface.html 182 import numba.cuda as nbcuda 183 desc = { 184 'shape': (length_pointer[0]), 185 'typestr': '>f8', 186 'data': (int(ffi.cast("intptr_t", array_pointer[0])), False), 187 'version': 2 188 } 189 # return read only Numba array 190 return nbcuda.from_cuda_array_interface(desc) 191 192 # Restore the Vector's data array 193 def restore_array(self): 194 """Restore an array obtained using get_array().""" 195 196 # Setup the pointer's pointer 197 array_pointer = ffi.new("CeedScalar **") 198 199 # libCEED call 200 err_code = lib.CeedVectorRestoreArray(self._pointer[0], array_pointer) 201 self._ceed._check_error(err_code) 202 203 # Restore an array obtained using getArrayRead 204 def restore_array_read(self): 205 """Restore an array obtained using get_array_read().""" 206 207 # Setup the pointer's pointer 208 array_pointer = ffi.new("CeedScalar **") 209 210 # libCEED call 211 err_code = lib.CeedVectorRestoreArrayRead( 212 self._pointer[0], array_pointer) 213 self._ceed._check_error(err_code) 214 215 @contextlib.contextmanager 216 def array(self, *shape, memtype=MEM_HOST): 217 """Context manager for array access. 218 219 Args: 220 shape (tuple): shape of returned numpy.array 221 **memtype: memory type of the array being passed, default CEED_MEM_HOST 222 223 224 Returns: 225 np.array: writable view of vector 226 227 Examples: 228 Constructing the identity inside a libceed.Vector: 229 230 >>> vec = ceed.Vector(16) 231 >>> with vec.array(4, 4) as x: 232 >>> x[...] = np.eye(4) 233 """ 234 x = self.get_array(memtype=memtype) 235 if shape: 236 x = x.reshape(shape) 237 yield x 238 self.restore_array() 239 240 @contextlib.contextmanager 241 def array_read(self, *shape, memtype=MEM_HOST): 242 """Context manager for read-only array access. 243 244 Args: 245 shape (tuple): shape of returned numpy.array 246 **memtype: memory type of the array being passed, default CEED_MEM_HOST 247 248 Returns: 249 np.array: read-only view of vector 250 251 Examples: 252 Constructing the identity inside a libceed.Vector: 253 254 >>> vec = ceed.Vector(6) 255 >>> vec.set_value(1.3) 256 >>> with vec.array_read(2, 3) as x: 257 >>> print(x) 258 """ 259 x = self.get_array_read(memtype=memtype) 260 if shape: 261 x = x.reshape(shape) 262 yield x 263 self.restore_array_read() 264 265 # Get the length of a Vector 266 def get_length(self): 267 """Get the length of a Vector. 268 269 Returns: 270 length: length of the Vector""" 271 272 length_pointer = ffi.new("CeedInt *") 273 274 # libCEED call 275 err_code = lib.CeedVectorGetLength(self._pointer[0], length_pointer) 276 self._ceed._check_error(err_code) 277 278 return length_pointer[0] 279 280 # Get the length of a Vector 281 def __len__(self): 282 """Get the length of a Vector. 283 284 Returns: 285 length: length of the Vector""" 286 287 length_pointer = ffi.new("CeedInt *") 288 289 # libCEED call 290 err_code = lib.CeedVectorGetLength(self._pointer[0], length_pointer) 291 self._ceed._check_error(err_code) 292 293 return length_pointer[0] 294 295 # Set the Vector to a given constant value 296 def set_value(self, value): 297 """Set the Vector to a constant value. 298 299 Args: 300 value: value to be used""" 301 302 # libCEED call 303 err_code = lib.CeedVectorSetValue(self._pointer[0], value) 304 self._ceed._check_error(err_code) 305 306 # Sync the Vector to a specified memtype 307 def sync_array(self, memtype=MEM_HOST): 308 """Sync the Vector to a specified memtype. 309 310 Args: 311 **memtype: memtype to be synced""" 312 313 # libCEED call 314 err_code = lib.CeedVectorSyncArray(self._pointer[0], memtype) 315 self._ceed._check_error(err_code) 316 317 # Compute the norm of a vector 318 def norm(self, normtype=NORM_2): 319 """Get the norm of a Vector. 320 321 Args: 322 **normtype: type of norm to be computed""" 323 324 norm_pointer = ffi.new("CeedScalar *") 325 326 # libCEED call 327 err_code = lib.CeedVectorNorm(self._pointer[0], normtype, norm_pointer) 328 self._ceed._check_error(err_code) 329 330 return norm_pointer[0] 331 332 # Take the reciprocal of a vector 333 def reciprocal(self): 334 """Take the reciprocal of a Vector.""" 335 336 # libCEED call 337 err_code = lib.CeedVectorReciprocal(self._pointer[0]) 338 self._ceed._check_error(err_code) 339 340 return self 341 342# ------------------------------------------------------------------------------ 343 344 345class _VectorWrap(Vector): 346 """Wrap a CeedVector pointer in a Vector object.""" 347 348 # Constructor 349 def __init__(self, ceed, pointer): 350 # CeedVector object 351 self._pointer = pointer 352 353 # Reference to Ceed 354 self._ceed = ceed 355 356# ------------------------------------------------------------------------------ 357