; mat_mat_64.asm  a direct coding of mat_mat.c 
;// mat_mat.c  compute matrix times matrix multiply
;#include <stdio.h>
;
;void mat_mat(long int n, double A[][n], double B[][], double C[][])
;{
;  long int i, j, k;
;
;  printf("mat_mat n=%ld \n",n);
;  for(i=0; i<n; i++)
;  {
;    for(j=0; j<n; j++)
;    {
;      C[i][j] = 0.0;
;      for(k=0; k<n; k++)
;      {
;        C[i][j] = C[i][j] + A[i][k]*B[k][j];
;      } // end k
;    } // end j
;  } // end i
;} // end mat_mat
;
;int main()
;{
;  long int n3 = 3;
;  double A1[3][3] = {{1.0,2.0,3.0},{4.0,5.0,6.0},{7.0,8.0,9.0}}; // A matrix
;  double B1[3][3] = {{2.0,3.0,4.0};{5.0,6.0,7.0},{8.0,9.0,1.0}}; // B matrix
;  double C1[3][3] = {{1.0,1.0,1.0};{1.0,1.0,1.0},{1.0,1.0,1.0}}; // C matrix;
;
;  printf("mat_mat_64.c running \n"); 
;  mat_mat(n3, A1, B1, C1);
;  printf("mat_mat returned =%f %f %f \n",C1[0][0], C1[0][1], C1[0][2]); 
;  printf("                  %f %f %f \n",C1[1][0], C1[1][1], C1[1][2]); 
;  printf("                  %f %f %f \n",C1[2][0], C1[2][1], C1[2][2]); 
;  return 0;
;} // end mat_mat.c

        extern printf

	section .data		; everything quad words, 8 bytes
dfmt3r:	db "mat_mat returned =%f %f %f", 10, 0 
fmtmsg:	db "mat_mat_64.asm running", 10, 0 
A1:	dq 1.0,2.0,3.0, 4.0,5.0,6.0, 7.0,8.0,9.0   ; A matrix
B1:	dq 2.0,3.0,4.0, 5.0,6.0,7.0, 8.0,9.0,1.0   ; B matrix
C1:	dq 1.0,1.0,1.0, 1.0,1.0,1.0, 1.0,1.0,1.0   ; C matrix;
x1:	dq 1.0,2.0,3.0          ; x vectro
n3:	dq 3			; length of vectors, size of matrix


	section .bss		; very unoptimized to get it correct
y1:	resq  3		        ; y vector result
	
	section .text
	global main
main:
        push    rbp		; set up stack frame, must be aligned

	mov	rdi,fmtmsg      ; address of format, standard register rdi
	mov	rax,0		; or can be  xor  rax,rax
        call    printf		; Call printf

	mov	rdi,[n3]	; set up to call mat_mat, value of n3
	mov	rsi,A1		; address of matrix A1
	mov	rdx,B1		; address of matrix B1
	mov	rcx,C1		; address of matrix C1
	call	mat_mat		; result returned in C1

	; print result matrix C1
	mov	rdi, dfmt3r
	movq	xmm0, qword[C1]
	movq	xmm1, qword[C1+8]
	movq	xmm2, qword[C1+16]
	mov	rax,3
	call	printf
	mov	rdi, dfmt3r
	movq	xmm0, qword[C1+24]
	movq	xmm1, qword[C1+32]
	movq	xmm2, qword[C1+40]
	mov	rax,3
	call	printf
	mov	rdi, dfmt3r
	movq	xmm0, qword[C1+48]
	movq	xmm1, qword[C1+56]
	movq	xmm2, qword[C1+64]
	mov	rax,3
	call	printf
	
	pop	rbp		; restore stack
	mov	rax,0		; normal, no error, return value
	ret			; return

; 	end main, start function
	
	section .bss
A:	resq 1			; address will be in rax
B:	resq 1			; address will be in rdx
C:	resq 1			; addres will be in rcx
i:	resq 1			; will be in r8  byte subscript [i]
j:	resq 1			; will be in r9  byte subscript [j]
k:	resq 1			; will be in r10  byte subscript [j]
n:	resq 1			; will be in r11  n
n8:	resq 1			; will be in r12 n*8
ij:	resq 1			; will be in r13 [i][j]
ik:	resq 1			; will be in r14 [i][k]
kj:	resq 1			; will be in r15 [k][j]
	
	section .data
zero:	dq 0.0
eight:	dq 8

	section .text
	global mat_mat		; mat_mat(n,A,B,C)
mat_mat:
        push	rbp		; save rbp

	mov	[n],rdi		; value of n
	mov	[A],rsi		; save address of A matrix
	mov	[B],rdx		; save address of B matrix
	mov	[C],rcx		; save address of C matric result

	mov	rax, [n]
	imul	qword [eight]
        mov	r12, rax	;  [*][]
	mov 	[n8], rax
	mov	r11, [n]	; n
	
	
	mov	r8, 0		; i = 0
iloop:				; for(i=0; i<n; i++)
	mov	[i],r8
	
	mov	r9, 0		; j = 0
jloop:				; for(j=0; j<n; j++)
	mov	r8, [i]
	mov	[j], r9

	mov	rax, [A]	; address of A matrix
	mov	rdx, [B]	; address of B matrix
	mov	rcx, [C]	; address of C matrix

	mov	rax, [i]	; i
	imul	qword [n]	; n
	add	rax, [j]	; j
	mov	r13,rax		; i*n+j  i*8*n+j*8 bytes

	fld	qword [zero]
	fstp	qword [rcx+r13]  ; C[i] = 0.0


	mov	r10, 0		; k = 0
kloop:				; for(k=0; k<n; k++)
	mov	r8, [i]
	mov	r9, [j]
	mov	[k], r10

	mov	rax, [i]	; i
	imul	qword [n]	; n
	add	rax, [k]	; k
	mov	r14,rax		; i*n+k  i*8*n+k*8 bytes

	mov	rax, [k]	; k
	imul	qword [n]	; n
	add	rax, [j]	; j
	mov	r15,rax		; k*n+j  i*8*n+j*8 bytes

	mov	rax, [A]	; address of A matrix
	mov	rdx, [B]	; address of x vector
	mov	rcx, [C]	; address of y vector

	fld     qword [rax+r14]     ; A[i][k]
	fmul	qword [rdx+r15]	    ; B[k][j]
	fadd	qword [rcx+r13]     ; C[i][j]
	fstp	qword [rcx+r13]	    ; C[i][j]

	mov	r10, [k]
	add	r10, 8
	mov	[k], r10
	cmp	r10, [n8]
	jl	kloop

	mov	r9, [j]
	add	r9, 8
	mov	[j], r9
	cmp	r9, [n8]
	jl	jloop

	mov	r8, [i]
	add	r8, 8
	mov	[i],r8
	cmp	r8, [n8]
	jl	iloop
	
        pop	rbp
        ret			; return
; end mat_mat_64.asm