FunctionLAPACKdsyrk
// 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 part of LAPACK: Linear Algebra Package, 
//  performs one of the symmetric rank k operations
// .
//    C := alpha*A*A**T + beta*C,   or    C := alpha*A**T*A + beta*C,
// .
// where alpha and beta are scalars, C is an n by n symmetric matrix
// and A is an n by k matrix in the first case and a k by n matrix
// in the second case.
// .
// reference:
//		https://netlib.org/lapack/explore-html/dc/d05/dsyrk_8f_source.html#l00168
library(title = 'FunctionLAPACKdsyrk')


// { 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 DSYRK
// 
//   =========== DOCUMENTATION ===========
// 
//  Online html documentation available at
//             http://www.netlib.org/lapack/explore-html/
// 
//   Definition:
//   ===========
// 
//        SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
// 
//        .. Scalar Arguments ..
//        DOUBLE PRECISION ALPHA,BETA
//        INTEGER K,LDA,LDC,N
//        CHARACTER TRANS,UPLO
//        ..
//        .. Array Arguments ..
//        DOUBLE PRECISION A(LDA,*),C(LDC,*)
//        ..
// 
// 
// \par Purpose:
//   =============
// >
// \verbatim
// >
// DSYRK  performs one of the symmetric rank k operations
// >
//    C := alpha*A*A**T + beta*C,
// >
// or
// >
//    C := alpha*A**T*A + beta*C,
// >
// where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
// and  A  is an  n by k  matrix in the first case and a  k by n  matrix
// in the second case.
// \endverbatim
// 
//   Arguments:
//   ==========
// 
// \param[in] UPLO
// \verbatim
//          UPLO is CHARACTER*1
//           On  entry,   UPLO  specifies  whether  the  upper  or  lower
//           triangular  part  of the  array  C  is to be  referenced  as
//           follows:
// >
//              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
//                                  is to be referenced.
// >
//              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
//                                  is to be referenced.
// \endverbatim
// >
// \param[in] TRANS
// \verbatim
//          TRANS is CHARACTER*1
//           On entry,  TRANS  specifies the operation to be performed as
//           follows:
// >
//              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
// >
//              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
// >
//              TRANS = 'C' or 'c'   C := alpha*A**T*A + beta*C.
// \endverbatim
// >
// \param[in] N
// \verbatim
//          N is INTEGER
//           On entry,  N specifies the order of the matrix C.  N must be
//           at least zero.
// \endverbatim
// >
// \param[in] K
// \verbatim
//          K is INTEGER
//           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
//           of  columns   of  the   matrix   A,   and  on   entry   with
//           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
//           of rows of the matrix  A.  K must be at least zero.
// \endverbatim
// >
// \param[in] ALPHA
// \verbatim
//          ALPHA is DOUBLE PRECISION.
//           On entry, ALPHA specifies the scalar alpha.
// \endverbatim
// >
// \param[in] A
// \verbatim
//          A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
//           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
//           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
//           part of the array  A  must contain the matrix  A,  otherwise
//           the leading  k by n  part of the array  A  must contain  the
//           matrix A.
// \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  TRANS = 'N' or 'n'
//           then  LDA must be at least  math.max( 1, n ), otherwise  LDA must
//           be at least  math.max( 1, k ).
// \endverbatim
// >
// \param[in] BETA
// \verbatim
//          BETA is DOUBLE PRECISION.
//           On entry, BETA specifies the scalar beta.
// \endverbatim
// >
// \param[in,out] C
// \verbatim
//          C is DOUBLE PRECISION array, dimension ( LDC, N )
//           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
//           upper triangular part of the array C must contain the upper
//           triangular part  of the  symmetric matrix  and the strictly
//           lower triangular part of C is not referenced.  On exit, the
//           upper triangular part of the array  C is overwritten by the
//           upper triangular part of the updated matrix.
//           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
//           lower triangular part of the array C must contain the lower
//           triangular part  of the  symmetric matrix  and the strictly
//           upper triangular part of C is not referenced.  On exit, the
//           lower triangular part of the array  C is overwritten by the
//           lower triangular part of the updated matrix.
// \endverbatim
// >
// \param[in] LDC
// \verbatim
//          LDC is INTEGER
//           On entry, LDC specifies the first dimension of C as declared
//           in  the  calling  (sub)  program.   LDC  must  be  at  least
//           math.max( 1, n ).
// \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
// >
//   =====================================================================
// end brief }

// @function performs one of the symmetric rank k operations
// .
//    C := alpha*A*A**T + beta*C,   or    C := alpha*A**T*A + beta*C,
// .
// where alpha and beta are scalars, C is an n by n symmetric matrix
// and A is an n by k matrix in the first case and a k by n matrix
// in the second case.
// .
// @param uplo  string        specifies whether the upper or lower triangular part of 
//                            the array C  is to be  referenced  as follows:
//              UPLO = 'U' or 'u'   Only the upper triangular part of  C is to be referenced.
//              UPLO = 'L' or 'l'   Only the lower triangular part of C is to be referenced.
//.
// @param trans string        specifies the operation to be performed as follows:
//              TRANS = 'N' or 'n'   C := alpha*A*A**T + beta*C.
//              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.
//              TRANS = 'C' or 'c'   C := alpha*A**T*A + beta*C.
// .
// @param n     int           specifies the order of the matrix C. N must be at least zero.
// @param k     int           On entry with:
//              TRANS = 'N' or 'n', K specifies the number of  columns   of  the   matrix   A.
//              TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number of rows of the matrix  A.  
//              K must be at least zero.
// .
// @param alpha float         scalar.
// @param a     matrix<float> matrix A.
// @param lda   int           specifies the first dimension of A.
// @param beta  float         scalar.
// @param c     matrix<float> matrix C, is overwritten by the lower triangular part of the updated matrix.
// @param ldc   int           specifies the first dimension of C
// @returns void, C is overwritten by the lower triangular part of the updated matrix.
export dsyrk (
  string uplo , string trans, 
  int    n    , int    k, 
  float  alpha, matrix<float> a    , int lda, 
  float  beta , matrix<float> c    , int ldc
  ) => // {
	int _n = n, int _n1 = _n - 1
	int _k = math.max(1, k)
	// Test the input parameters.
	int  _nrowa = lsame(trans,'N') ? n : _k
	bool _upper = lsame(uplo,'U')
	// 
	int _info = 0
	switch
		(not _upper) and (not lsame(uplo, 'L')) 									=> _info := 1
		(not lsame(trans,'N')) and (not lsame(trans,'T')) and (not lsame(trans,'C'))=> _info := 2
		_n < 0 																		=> _info := 3
		k < 0 																		=> _info := 4
		lda < math.max(1, _nrowa) 													=> _info := 7
		ldc < math.max(1, _n) 														=> _info := 10
	if (_info != 0)
		xerbla('DSYRK ', _info) // ### log and break, TODO: can build in messages directly..
	switch
		// Quick return if possible.
		(_n == 0 or ((alpha == 0.0 or k == 0) and beta == 1.0)) => 0 // RETURN
		alpha == 0.0 => 
			if _upper
				if beta == 0.0
					0
					for _j = 0 to _n1
						for _i = 0 to _j
							matrix.set(c, _i, _j, 0.0)
				else
					for _j = 0 to _n1
						for _i = 0 to _j
							matrix.set(c, _i, _j, beta * matrix.get(c, _i, _j))
        	else
				if beta == 0.0
					for _j = 0 to _n1
						for _i = _j to _n1
							matrix.set(c, _i, _j, 0.0)
				else
					for _j = 0 to _n1
						for _i = _j to _n1
							matrix.set(c, _i, _j, beta * matrix.get(c, _i, _j))
        	0 // RETURN
		=>
			// Start the operations.
			if lsame(trans, 'N')
				// Form  C := alpha*A*A**T + beta*C.
				if _upper
					for _j = 0 to _n1
						if beta == 0.0
							for _i = 0 to _j
								matrix.set(c, _i, _j, 0.0)
						else if beta != 1.0
							for _i = 0 to _j
								matrix.set(c, _i, _j, beta * matrix.get(c, _i, _j))
						for _l = 0 to _k-1 // ### modified logic on _k
							if matrix.get(a, _j, _l) != 0.0
								float _temp = alpha * matrix.get(a, _j, _l)
								for _i = 0 to _j
									matrix.set(c, _i, _j, matrix.get(c, _i, _j) + _temp * matrix.get(a, _i, _l))
				else
					for _j = 0 to _n1
						if beta == 0.0
							for _i = _j to _n1
								matrix.set(c, _i, _j, 0.0)
						else if beta != 1.0
							for _i = _j to _n1
								matrix.set(c, _i, _j, beta * matrix.get(c, _i, _j))
						for _l = 0 to _k-1 // ### modified logic on _k
							if matrix.get(a, _j, _l) != 0.0
								float _temp = alpha * matrix.get(a, _j, _l)
								for _i = _j to _n1
									matrix.set(c, _i, _j, matrix.get(c, _i, _j) + _temp * matrix.get(a, _i, _l))
			else
				// Form  C := alpha*A**T*A + beta*C.
				if _upper
					for _j = 0 to _n1
						for _i = 0 to _j
							float _temp = 0.0
							for _l = 0 to _k-1 // ### modified logic on _k
								_temp += matrix.get(a, _l, _i) * matrix.get(a, _l, _j)
							if beta == 0.0
								matrix.set(c, _i, _j, alpha * _temp)
							else
								matrix.set(c, _i, _j, alpha * _temp + beta * matrix.get(c, _i, _j))
				else
					for _j = 0 to _n1
						for _i = _j to _n1
							float _temp = 0.0
							for _l = 0 to _k-1 // ### modified logic on _k
								_temp += matrix.get(a, _l, _i) * matrix.get(a, _l, _j)
							if beta == 0.0
								matrix.set(c, _i, _j, alpha * _temp)
							else
								matrix.set(c, _i, _j, alpha * _temp + beta * matrix.get(c, _i, _j))
			0 // RETURN
	// End of DSYRK
// { example:
// reference:
//		https://www.ibm.com/docs/en/essl/6.1?topic=mos-ssyrk-dsyrk-csyrk-zsyrk-cherk-zherk-rank-update-real-complex-symmetric-complex-hermitian-matrix
A = matrix.new<float>(8, 0)
matrix.add_col(A, 0, array.from(0.0, 1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0))
matrix.add_col(A, 1, array.from(8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0))
C = matrix.new<float>(0, 8)
matrix.add_row(C, 0, array.from(      0.0,       1.0,       3.0,       6.0,      10.0,      15.0,      21.0, 28.0))
matrix.add_row(C, 1, array.from(float(na),       2.0,       4.0,       7.0,      11.0,      16.0,      22.0, 29.0))
matrix.add_row(C, 2, array.from(float(na), float(na),       5.0,       8.0,      12.0,      17.0,      23.0, 30.0))
matrix.add_row(C, 3, array.from(float(na), float(na), float(na),       9.0,      13.0,      18.0,      24.0, 31.0))
matrix.add_row(C, 4, array.from(float(na), float(na), float(na), float(na),      14.0,      19.0,      25.0, 32.0))
matrix.add_row(C, 5, array.from(float(na), float(na), float(na), float(na), float(na),      20.0,      26.0, 33.0))
matrix.add_row(C, 6, array.from(float(na), float(na), float(na), float(na), float(na), float(na),      27.0, 34.0))
matrix.add_row(C, 7, array.from(float(na), float(na), float(na), float(na), float(na), float(na), float(na), 35.0))

import RicardoSantos/DebugConsole/10 as console
[__T, __C] = console.init()
console.queue_one(__C, 'A:\n' + str.tostring(A))
console.queue_one(__C, 'C:\n' + str.tostring(C))
dsyrk('U', 'N', 8, 2,  1.0, A, 9, 1.0, C, 10)
console.queue_one(__C, 'solution X:\n' + str.tostring(C))

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

