Merge branch 'master' into script_devel to update with alternative inversion methods (ebfa84bb) · Commits · Giacomo Mulas / NP_TMcode

build/configure.sh

+27 −1

Original line number	Diff line number	Diff line
		@@ -16,6 +16,7 @@ LAPACK="auto"
		LDFLAGS=""
		LIBMODE="static"
		MAGMA="auto"
		MAGMA_INVERT_FLAGS=""
		NVTXFLAGS=""
		OMPMODE="auto"
		OFFLOAD="auto"
		@@ -98,6 +99,7 @@ function print_help {
		echo " than once). "
		echo "--with-magma=[MAGMA] Use specified MAGMA distribution (DEFAULT). "
		echo "--without-magma Disable MAGMA. "
		echo "--enable-magma-invert=OPTION Use specified MAGMA matrix inversion function (options = auto\|getri\|gesv\|rbt, default getri). "
		echo " "
		echo "Some influential environment variables are: "
		echo " "
		@@ -162,7 +164,7 @@ do
		if [ "x$DEBUGFLAGS" = "x" ]; then
		DEBUGFLAGS=" -DDEBUG_${dbg_feature}"
		else
		DEBUGFLAGS="${DEBUGFLAGS} -DDEBUG_{dbg_feature}"
		DEBUGFLAGS="${DEBUGFLAGS} -DDEBUG_${dbg_feature}"
		fi
		fi
		elif [ "x$cut_arg" = "x--enable-fortran" ]; then
		@@ -220,6 +222,18 @@ do
		MAGMA_HOME=$(echo $arg \| cut -d '=' -f2)
		elif [ "x$cut_arg" = "x--without-magma" ]; then
		MAGMA="no"
		elif [ "x$cut_arg" = "x--enable-magma-invert" ]; then
		MAGMA_INVERT_CHOICE=$(echo $arg \| cut -d '=' -f2)
		if [ "x$MAGMA_INVERT_CHOICE" = "xgetri" -o "x$MAGMA_INVERT_CHOICE" = "xauto" ]; then
		MAGMA_INVERT_FLAGS=""
		elif [ "x$MAGMA_INVERT_CHOICE" = "xgesv" ]; then
		MAGMA_INVERT_FLAGS=" -DUSE_ZGESV_GPU"
		elif [ "x$MAGMA_INVERT_CHOICE" = "xrbt" ]; then
		MAGMA_INVERT_FLAGS=" -DUSE_ZGESV_RBT"
		else
		echo "ERROR: unrecognized --enable-magma-invert option \"$MAGMA_INVERT_CHOICE\""
		exit 1
		fi
		else
		echo "ERROR: unrecognized argument \"$arg\""
		exit 1
		@@ -732,6 +746,11 @@ if [ "x$MAGMA" != "xno" ]; then
		exit 2
		fi
		else
		if [ "x$MAGMAFLAGS" = "x" ]; then
		MAGMAFLAGS="$MAGMA_INVERT_FLAGS"
		else
		MAGMAFLAGS="$MAGMAFLAGS$MAGMA_INVERT_FLAGS"
		fi
		echo "yes"
		fi
		else
		@@ -832,6 +851,13 @@ if [ "x${MAGMAFLAGS}" = "x" ]; then
		echo "INFO: MAGMA is disabled."
		else
		echo "INFO: MAGMA is enabled."
		if [ "x${MAGMA_INVERT_FLAGS}" = "x" ]; then
		echo "INFO: using LU factorisation for matrix inversion."
		elif [ "x${MAGMA_INVERT_FLAGS}" = "x -DUSE_ZGESV_GPU" ]; then
		echo "INFO: using MAGMA zgesv_gpu function for matrix inversion."
		elif [ "x${MAGMA_INVERT_FLAGS}" = "x -DUSE_ZGESV_RBT" ]; then
		echo "INFO: using MAGMA zgesv_rbt function for matrix inversion."
		fi
		fi
		if [ "x${NVTXFLAGS}" = "x" ]; then
		echo "INFO: NVTX profiling is disabled."

src/cluster/cluster.cpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -841,7 +841,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
		#endif
		// we put the accuracygoal in, get the actual accuracy back out
		double actualaccuracy = cid->accuracygoal;
		invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, mxndm, cid->proc_device);
		invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);
		// in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not
		#ifdef USE_REFINEMENT
		if (cid->refinemode==2) {

src/include/algebraic.h

+3 −1

Original line number	Diff line number	Diff line
		@@ -29,6 +29,8 @@
		#ifndef INCLUDE_ALGEBRAIC_H_
		#define INCLUDE_ALGEBRAIC_H_

		using namespace std;

		/*! \brief Perform in-place matrix inversion.
		*
		* \param mat: `complex double **` The matrix to be inverted (must be a square matrix).
		@@ -38,6 +40,6 @@
		* optional, defaults to 0).
		* \param target_device: `int` ID of target GPU, if available (defaults to 0).
		*/
		void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, np_int max_size=0, int target_device=0);
		void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, const string& output_path, int jxi488, np_int max_size=0, int target_device=0);

		#endif

src/include/magma_calls.h

+30 −4

Original line number	Diff line number	Diff line
		@@ -20,13 +20,14 @@
		*
		*/

		#include <string>

		#ifndef INCLUDE_MAGMA_CALLS_H_
		#define INCLUDE_MAGMA_CALLS_H_

		/*! \brief Invert a complex matrix with double precision elements.
		*
		* Use MAGMA to perform an in-place matrix inversion for a complex
		* matrix with double precision elements.
		* call magma_zinvert1() to do the actual inversion
		*
		* \param mat: Matrix of complex. The matrix to be inverted.
		* \param n: `np_int` The number of rows and columns of the [n x n] matrix.
		@@ -35,11 +36,22 @@
		*/
		void magma_zinvert(dcomplex **mat, np_int n, int &jer, int device_id=0);

		/*! \brief Invert a complex matrix with double precision elements, applying iterative refinement of the solution
		/*! \brief Auxiliary function to Invert a complex matrix with double precision elements.
		*
		* Use MAGMA to perform an in-place matrix inversion for a complex
		* matrix with double precision elements.
		*
		* \param inva: reference to the pointer to the first element of the matrix to be inverted on entry, inverted matrix on exit.
		* \param n: `np_int` The number of rows and columns of the [n x n] matrix.
		* \param jer: `int &` Reference to an integer return flag.
		* \param device_id: `int` ID of the device for matrix inversion offloading.
		*/
		void magma_zinvert1(dcomplex * &inva, np_int n, int &jer, int device_id);

		/*! \brief Invert a complex matrix with double precision elements, applying iterative refinement of the solution
		*
		* call magma_zinvert1() to perform the first matrix inversion, then magma_refine() to do the refinement (only if maxrefiters is >0)
		*
		* \param mat: Matrix of complex. The matrix to be inverted.
		* \param n: `np_int` The number of rows and columns of the [n x n] matrix.
		* \param jer: `int &` Reference to an integer return flag.
		@@ -47,6 +59,20 @@ void magma_zinvert(dcomplex **mat, np_int n, int &jer, int device_id=0);
		* \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy
		* \param device_id: `int` ID of the device for matrix inversion offloading.
		*/
		void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id);
		void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);

		/*! \brief apply iterative refinement of the solution of a matrix inversion
		*
		* iteratively compute and apply a correction to the inverse inva of the complex matrix aorig, for a maximum number of maxiters times, or until achieving a maximum residual better than accuracygoal
		*
		* \param aorig: pointer to the first element of the matrix of complex to be inverted.
		* \param inva: pointer to the first element of inverse.
		* \param n: `np_int` The number of rows and columns of the [n x n] matrices.
		* \param jer: `int &` Reference to an integer return flag.
		* \param maxrefiters: `int` Maximum number of refinement iterations to apply.
		* \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy
		* \param device_id: `int` ID of the device for matrix inversion offloading.
		*/
		void magma_refine(dcomplex aorig, dcomplex inva, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);

		#endif

src/inclusion/inclusion.cpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -1480,7 +1480,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
		#endif
		// we the accuracygoal in, get the actual accuracy back out
		double actualaccuracy = cid->accuracygoal;
		invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, mxndm, cid->proc_device);
		invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);
		// in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not
		#ifdef USE_REFINEMENT
		if (cid->refinemode==2) {