add support for matrix inversion via Cholesky decomposition

CMakeLists.txt

@@ -1,8 +1,8 @@
 # Translate BLAS/LAPACK subset to C++. Compile and setup tests.
 
-cmake_minimum_required(VERSION 3.16..3.25)
+cmake_minimum_required(VERSION 3.16..3.30)
 
-project(linalg VERSION 1.1
+project(linalg VERSION 1.2
   DESCRIPTION "BLAS/LAPACK functions translated to C++ as required by the LAMMPS MD package"
   LANGUAGES CXX)
 

fortran/dlauu2.f

@@ -0,0 +1,198 @@
+*> \brief \b DLAUU2 computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm).
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLAUU2 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlauu2.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlauu2.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlauu2.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DLAUU2( UPLO, N, A, LDA, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   A( LDA, * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DLAUU2 computes the product U * U**T or L**T * L, where the triangular
+*> factor U or L is stored in the upper or lower triangular part of
+*> the array A.
+*>
+*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
+*> overwriting the factor U in A.
+*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
+*> overwriting the factor L in A.
+*>
+*> This is the unblocked form of the algorithm, calling Level 2 BLAS.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the triangular factor stored in the array A
+*>          is upper or lower triangular:
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the triangular factor U or L.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension (LDA,N)
+*>          On entry, the triangular factor U or L.
+*>          On exit, if UPLO = 'U', the upper triangle of A is
+*>          overwritten with the upper triangle of the product U * U**T;
+*>          if UPLO = 'L', the lower triangle of A is overwritten with
+*>          the lower triangle of the product L**T * L.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0: successful exit
+*>          < 0: if INFO = -k, the k-th argument had an illegal value
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date December 2016
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE DLAUU2( UPLO, N, A, LDA, INFO )
+*
+*  -- LAPACK auxiliary routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     December 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I
+      DOUBLE PRECISION   AII
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   DDOT
+      EXTERNAL           LSAME, DDOT
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEMV, DSCAL, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DLAUU2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Compute the product U * U**T.
+*
+         DO 10 I = 1, N
+            AII = A( I, I )
+            IF( I.LT.N ) THEN
+               A( I, I ) = DDOT( N-I+1, A( I, I ), LDA, A( I, I ), LDA )
+               CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
+     $                     LDA, A( I, I+1 ), LDA, AII, A( 1, I ), 1 )
+            ELSE
+               CALL DSCAL( I, AII, A( 1, I ), 1 )
+            END IF
+   10    CONTINUE
+*
+      ELSE
+*
+*        Compute the product L**T * L.
+*
+         DO 20 I = 1, N
+            AII = A( I, I )
+            IF( I.LT.N ) THEN
+               A( I, I ) = DDOT( N-I+1, A( I, I ), 1, A( I, I ), 1 )
+               CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
+     $                     A( I+1, I ), 1, AII, A( I, 1 ), LDA )
+            ELSE
+               CALL DSCAL( I, AII, A( I, 1 ), LDA )
+            END IF
+   20    CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of DLAUU2
+*
+      END