FunctionLAPACKdtrsm
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5

// @description subroutine in the LAPACK:linear algebra package, used to solve one of the following matrix equations:
// 
//     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
// 
//  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
//  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
// 
//     op( A ) = A   or   op( A ) = A**T.
// 
//  The matrix X is overwritten on B.
//
// reference:
//		https://netlib.org/lapack/explore-html/de/da7/dtrsm_8f_source.html#l00180
library(title='FunctionLAPACKdtrsm')

// { helper functions:

// https://netlib.org/lapack/explore-html/d0/dbf/_b_l_a_s_2_s_r_c_2lsame_8f_source.html#l00052
// test if letter in 'a' equals 'b' regardless of case.
lsame (string a, string b) =>
    str.lower(a) == str.lower(b)
//

// https://netlib.org/lapack/explore-html/d1/dc0/_b_l_a_s_2_s_r_c_2xerbla_8f_source.html#l00059
xerbla (string srname, int info) =>
    runtime.error(str.format(' ** On entry to {0} parameter number {1} had an illegal value', srname, info))

// end: helper functions }


// { brief: DTRSM
// 
// =========== DOCUMENTATION ===========
// 
//  Online html documentation available at
//           http://www.netlib.org/lapack/explore-html/
// 
// Definition:
// ===========
// 
//      SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
// 
//      .. Scalar Arguments ..
//      DOUBLE PRECISION ALPHA
//      INTEGER LDA,LDB,M,N
//      CHARACTER DIAG,SIDE,TRANSA,UPLO
//      ..
//      .. Array Arguments ..
//      DOUBLE PRECISION A(LDA,*),B(LDB,*)
//      ..
// 
// 
//  \par Purpose:
// =============
// 
//  \verbatim
// 
//  DTRSM  solves one of the matrix equations
// 
//     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
// 
//  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
//  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
// 
//     op( A ) = A   or   op( A ) = A**T.
// 
//  The matrix X is overwritten on B.
//  \endverbatim
// 
// Arguments:
// ==========
// 
//  \param[in] SIDE
//  \verbatim
//           SIDE is CHARACTER*1
//            On entry, SIDE specifies whether op( A ) appears on the left
//            or right of X as follows:
// 
//               SIDE = 'L' or 'l'   op( A )*X = alpha*B.
// 
//               SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
//  \endverbatim
// 
//  \param[in] UPLO
//  \verbatim
//           UPLO is CHARACTER*1
//            On entry, UPLO specifies whether the matrix A is an upper or
//            lower triangular matrix as follows:
// 
//               UPLO = 'U' or 'u'   A is an upper triangular matrix.
// 
//               UPLO = 'L' or 'l'   A is a lower triangular matrix.
//  \endverbatim
// 
//  \param[in] TRANSA
//  \verbatim
//           TRANSA is CHARACTER*1
//            On entry, TRANSA specifies the form of op( A ) to be used in
//            the matrix multiplication as follows:
// 
//               TRANSA = 'N' or 'n'   op( A ) = A.
// 
//               TRANSA = 'T' or 't'   op( A ) = A**T.
// 
//               TRANSA = 'C' or 'c'   op( A ) = A**T.
//  \endverbatim
// 
//  \param[in] DIAG
//  \verbatim
//           DIAG is CHARACTER*1
//            On entry, DIAG specifies whether or not A is unit triangular
//            as follows:
// 
//               DIAG = 'U' or 'u'   A is assumed to be unit triangular.
// 
//               DIAG = 'N' or 'n'   A is not assumed to be unit
//                                   triangular.
//  \endverbatim
// 
//  \param[in] M
//  \verbatim
//           M is INTEGER
//            On entry, M specifies the number of rows of B. M must be at
//            least zero.
//  \endverbatim
// 
//  \param[in] N
//  \verbatim
//           N is INTEGER
//            On entry, N specifies the number of columns of B.  N must be
//            at least zero.
//  \endverbatim
// 
//  \param[in] ALPHA
//  \verbatim
//           ALPHA is DOUBLE PRECISION.
//            On entry,  ALPHA specifies the scalar  alpha. When  alpha is
//            zero then  A is not referenced and  B need not be set before
//            entry.
//  \endverbatim
// 
//  \param[in] A
//  \verbatim
//           A is DOUBLE PRECISION array, dimension ( LDA, k ),
//            where k is m when SIDE = 'L' or 'l'
//              and k is n when SIDE = 'R' or 'r'.
//            Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
//            upper triangular part of the array  A must contain the upper
//            triangular matrix  and the strictly lower triangular part of
//            A is not referenced.
//            Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
//            lower triangular part of the array  A must contain the lower
//            triangular matrix  and the strictly upper triangular part of
//            A is not referenced.
//            Note that when  DIAG = 'U' or 'u',  the diagonal elements of
//            A  are not referenced either,  but are assumed to be  unity.
//  \endverbatim
// 
//  \param[in] LDA
//  \verbatim
//           LDA is INTEGER
//            On entry, LDA specifies the first dimension of A as declared
//            in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
//            LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
//            then LDA must be at least max( 1, n ).
//  \endverbatim
// 
//  \param[in,out] B
//  \verbatim
//           B is DOUBLE PRECISION array, dimension ( LDB, N )
//            Before entry,  the leading  m by n part of the array  B must
//            contain  the  right-hand  side  matrix  B,  and  on exit  is
//            overwritten by the solution matrix  X.
//  \endverbatim
// 
//  \param[in] LDB
//  \verbatim
//           LDB is INTEGER
//            On entry, LDB specifies the first dimension of B as declared
//            in  the  calling  (sub)  program.   LDB  must  be  at  least
//            max( 1, m ).
//  \endverbatim
// 
// Authors:
// ========
// 
//  \author Univ. of Tennessee
//  \author Univ. of California Berkeley
//  \author Univ. of Colorado Denver
//  \author NAG Ltd.
// 
//  \ingroup double_blas_level3
// 
//  \par Further Details:
// =====================
// 
//  \verbatim
// 
//   Level 3 Blas routine.
// 
// 
//   -- Written on 8-February-1989.
//      Jack Dongarra, Argonne National Laboratory.
//      Iain Duff, AERE Harwell.
//      Jeremy Du Croz, Numerical Algorithms Group Ltd.
//      Sven Hammarling, Numerical Algorithms Group Ltd.
//  \endverbatim
// 
// =====================================================================
// }

// @function solves one of the matrix equations
//     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
//  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
//  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
//     op( A ) = A   or   op( A ) = A**T.
//  The matrix X is overwritten on B.
//
// @param side   string       , On entry, SIDE specifies whether op( A ) appears on the left or right of X as follows:
//               SIDE = 'L' or 'l'   op( A )*X = alpha*B.
//               SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
// @param uplo   string       , specifies whether the matrix A is an upper or lower triangular matrix as follows:
//               UPLO = 'U' or 'u'   A is an upper triangular matrix.
//               UPLO = 'L' or 'l'   A is a lower triangular matrix.
// @param transa string       , specifies the form of op( A ) to be used in the matrix multiplication as follows:
//               TRANSA = 'N' or 'n'   op( A ) = A.
//               TRANSA = 'T' or 't'   op( A ) = A**T.
//               TRANSA = 'C' or 'c'   op( A ) = A**T.
// @param diag   string       , specifies whether or not A is unit triangular as follows:
//               DIAG = 'U' or 'u'   A is assumed to be unit triangular.
//               DIAG = 'N' or 'n'   A is not assumed to be unit triangular.
// @param m      int          , the number of rows of B. M must be at least zero.
// @param n      int          , the number of columns of B. N must be at least zero.
// @param alpha  float        , specifies the scalar alpha. When alpha is zero then A is not referenced and B need not be set before entry.
// @param a      matrix<float>, Triangular matrix.
// @param lda    int          , specifies the first dimension of A.
// @param b      matrix<float>, right-hand side matrix B, and on exit is overwritten by the solution matrix X.
// @param ldb    int          , specifies the first dimension of B.
// @returns void, modifies matrix b.
// usage:
//		dtrsm ('L', 'U', 'N', 'N', 5, 3, 1.0, a, 7, b, 6)
export dtrsm (
  string side, string uplo, string transa, string diag,
  int m, int n, float alpha, 
  matrix<float> a, int lda, 
  matrix<float> b, int ldb
  ) => // {
	int _m = m, int _m1 = _m - 1
	int _n = n, int _n1 = _n - 1
	// Test the input parameters.
    bool _lside  = lsame(side, 'L')
    int  _nrowa  = _lside ? _m : _n
    bool _nounit = lsame(diag, 'N')
    bool _upper  = lsame(uplo, 'U')
	//
	int _info = 0
	switch
		(not _lside) and (not lsame(side,'R')) 											=> _info := 1
		(not _upper) and (not lsame(uplo,'L')) 											=> _info := 2
		(not lsame(transa,'N')) and (not lsame(transa,'T')) and (not lsame(transa,'C')) => _info := 3
		(not lsame(diag,'U')) and (not lsame(diag,'N')) 								=> _info := 4
		_m < 0 																			=> _info := 5
		_n < 0 																			=> _info := 6
		lda < math.max(1, _nrowa) 														=> _info := 9
		ldb < math.max(1, _m) 															=> _info := 11
	if _info != 0
		xerbla('DTRSM ', _info)// ### log and break, TODO: can build in messages directly..
	switch
		// Quick return if possible.
		_m == 0 or _n == 0 	=> 0 // RETURN
		alpha == 0.0		=>
			for _j = 0 to _n1
				for _i = 0 to _m1
					matrix.set(b, _i, _j, 0.0)
			0 // RETURN
		=>
			//    Start the operations.
			if _lside
				if lsame(transa, 'N')
					// Form  B := alpha*inv( A )*B.
					if _upper
						for _j = 0 to _n1
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _j, alpha * matrix.get(b, _i, _j))
							for _k = _m1 to 0
								float _Bkj = matrix.get(b, _k, _j)
								if _Bkj != 0.0
									if _nounit
										matrix.set(b, _k, _j, _Bkj / matrix.get(a, _k, _k))
									if _k > 0 // ### protection against index out of bounds.
										for _i = 0 to _k - 1 // ### potencial out of bounds.
										//DO 40 i = 1,k - 1
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - matrix.get(b, _k, _j) * matrix.get(a, _i, _k))
					else
						for _j = 0 to _n1
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _j, alpha * matrix.get(b, _i, _j))
							for _k = 0 to _m1
								if matrix.get(b, _k, _j) != 0.0
									if _nounit
										matrix.set(b, _k, _j, matrix.get(b, _k, _j) / matrix.get(a, _k, _k))
									if _k > 0 // ### protection against index out of bounds.
										for _i = _k + 1 to _m1 // ### potencial out of bounds.
										//DO 80 i = k + 1,m
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - matrix.get(b, _k, _j) * matrix.get(a, _i, _k))
				else
					// Form  B := alpha*inv( A**T )*B.
					if _upper
						for _j = 0 to _n1
							for _i = 0 to _m1
								float _temp = alpha * matrix.get(b, _i, _j)
								if _i > 0 // ### protection against index out of bounds.
									for _k = 0 to _i - 1 // ### potential out of bounds.
									//DO 110 k = 1,i - 1
										_temp -= matrix.get(a, _k, _i) * matrix.get(b, _k, _j)
								if _nounit
									_temp /= matrix.get(a, _i, _i)
								matrix.set(b, _i, _j, _temp)
					else
						for _j = 0 to _n1
							for _i = _m1 to 0
								float _temp = alpha * matrix.get(b, _i, _j)
								if _i < _m1 // ### protection against index out of bounds.
									for _k = _i + 1 to _m1 // ### potential out of bounds.
									//DO 140 k = i + 1,m
										_temp -= matrix.get(a, _k, _i) * matrix.get(b, _k, _j)
								if _nounit
									_temp /= matrix.get(a, _i, _i)
								matrix.set(b, _i, _j, _temp)
			else
				if lsame(transa, 'N')
					// Form  B := alpha*B*inv( A ).
					if _upper
						for _j = 0 to _n1
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _j, alpha * matrix.get(b, _i, _j))
							if _j > 0 // ### protection against index out of bounds.
								for _k = 0 to _j - 1 // ### potential out of bounds.
								//DO 190 k = 1,j - 1
									if matrix.get(a, _k, _j) != 0.0
										for _i = 0 to _m1
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - matrix.get(a, _k, _j) * matrix.get(b, _i, _k))
							if _nounit
								float _temp = 1.0 / matrix.get(a, _j, _j)
								for _i = 0 to _m1
									matrix.set(b, _i, _j, _temp * matrix.get(b, _i, _j))
					else
						for _j = _n1 to 0
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _j, alpha * matrix.get(b, _i, _j))
							if _j < _n1 // ### protection against index out of bounds.
								for _k = _j + 1 to _n1 // ### potential out of bounds.
								//DO 240 k = j + 1,n
									if matrix.get(a, _k, _j) != 0.0
										for _i = 0 to _m1
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - matrix.get(a, _k, _j) * matrix.get(b, _i, _k))
							if _nounit
								float _temp = 1.0 / matrix.get(a, _j, _j)
								for _i = 0 to _m1
									matrix.set(b, _i, _j, _temp * matrix.get(b, _i, _j))
				else
					// Form  B := alpha*B*inv( A**T ).
					if _upper
						for _k = _n1 to 0
							if _nounit
								float _temp = 1.0 / matrix.get(a, _k, _k)
								for _i = 0 to _m1
									matrix.set(b, _i, _k, _temp * matrix.get(b, _i, _k))
							if _k > 0 // ### protection against index out of bounds.
								for _j = 0 to _k - 1 // ### potential out of bounds.
								//DO 290 j = 1,k - 1
									if matrix.get(a, _j, _k) != 0.0
										float _temp = matrix.get(a, _j, _k)
										for _i = 0 to _m1
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - _temp * matrix.get(b, _i, _k))
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _k, alpha * matrix.get(b, _i, _k))
					else
						for _k = 0 to _n1
							if _nounit
								float _temp = 1.0 / matrix.get(a, _k, _k)
								for _i = 0 to _m1
									matrix.set(b, _i, _k, _temp * matrix.get(b, _i, _k))
							if _k < _n1 // ### protection against index out of bounds.
								for _j = _k + 1 to _n1 // ### potential out of bounds.
								//DO 340 j = k + 1,n
									if matrix.get(a, _j, _k) != 0.0
										float _temp = matrix.get(a, _j, _k)
										for _i = 0 to _m1
											matrix.set(b, _i, _j, matrix.get(b, _i, _j) - _temp * matrix.get(b, _i, _k))
							if alpha != 1.0
								for _i = 0 to _m1
									matrix.set(b, _i, _k, alpha * matrix.get(b, _i, _k))
			0 // RETURN
	// End of DTRSM
// { example:
// tested only on this example..

// string side = input.string(defval = 'L', options = ['L', 'l', 'R', 'r'])
// string uplo = input.string(defval = 'U', options = ['U', 'u', 'L', 'l'])
// string transa = input.string(defval = 'N', options = ['N', 'n', 'T', 't', 'C', 'c'])
// string diag = input.string(defval = 'N', options = ['U', 'u', 'N', 'n'])
// int m = input.int(5)
// int n = input.int(3)
// float alpha = input.float(1.0)
// int k = switch side
// 	'L' => m
// 	'R' => n
// lda = switch side
// 	'L' => m
// 	'R' => n
// ldb = m


a = matrix.new<float>(0, 5)//(lda, k, 1.0)
matrix.add_row(a, 0, array.from(3.0, -1.0, 2.0, 2.0, 1.0))
matrix.add_row(a, 1, array.from(float(na), -2.0, 4.0, -1.0, 3.0))
matrix.add_row(a, 2, array.from(float(na), float(na), -3.0, 0.0, 2.0))
matrix.add_row(a, 3, array.from(float(na), float(na), float(na), 4.0, -2.0))
matrix.add_row(a, 4, array.from(float(na), float(na), float(na), float(na), 1.0))
// matrix.add_row(a, 5), matrix.add_row(a, 6)
b = matrix.new<float>(0, 3)//(ldb, n, 1.0)
matrix.add_row(b, 0, array.from(6.0, 10.0, -2.0))
matrix.add_row(b, 1, array.from(-16.0, -1.0, 6.0))
matrix.add_row(b, 2, array.from(-2.0, 1.0, -4.0))
matrix.add_row(b, 3, array.from(14.0, 0.0, -14.0))
matrix.add_row(b, 4, array.from(-1.0, 2.0, 1.0))
// matrix.add_row(b, 5)
import RicardoSantos/DebugConsole/10 as console
[__T, __C] = console.init()
console.queue_one(__C, 'A:\n' + str.tostring(a))
console.queue_one(__C, 'B:\n' + str.tostring(b))
//dtrsm (side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
// https://www.ibm.com/docs/en/essl/6.1?topic=mos-strsm-dtrsm-ctrsm-ztrsm-solution-triangular-systems-equations-multiple-right-hand-sides
// example 1:
dtrsm ('L', 'U', 'N', 'N', 5, 3, 1.0, a, 7, b, 6)
console.queue_one(__C, 'solution X:\n' + str.tostring(b))

console.update(__T, __C)
//}
//}