; mat_vec_64.asm  a direct coding of mat_vec.c
;// mat_vec.c  compute matrix times vector multiply
;#include <stdio.h>
;
;void mat_vec(long int n, double A[][n], double x[], double y[])
;{
;  long int i, j;
;
;  printf("mat_vec n=%ld \n",n);
;  for(i=0; i<n; i++)
;  {
;    y[i] = 0.0;
;    for(j=0; j<n; j++)
;    {
;      y[i] = y[i] + A[i][j]*x[j];
;    } // end j
;  } // end i
;} // end mat_vec
;
;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 x1[3] = {1.0,2.0,3.0};          // x vector
;  double y1[3];
;
;  printf("mat_vec_64.asm running \n"); 
;  mat_vec(n3, A1, x1, y1);
;  printf("mat_vec returned =%f %f %f \n",y1[0], y1[1], y1[2]); 
;  return 0;
;} // end mat_vec.c

        extern printf

	section .data		; everything quad words, 8 bytes
dfmt2r:	db "mat_vec returned =%f %f %f", 10, 0 
fmtmsg:	db "mat_vec_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
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 C function

	mov	rdi,[n3]	; set up to call mat_vec, value of n3
	mov	rsi,A1		; address of matrix A1
	mov	rdx,x1		; address of vector x1
	mov	rcx,y1		; address of vector y1
	call	mat_vec		; result returned in y1

	; print result vector y1 
	mov	rdi, dfmt2r
	movq	xmm0, qword[y1]
	movq	xmm1, qword[y1+8]
	movq	xmm2, qword[y1+16]
	mov	rax,3
	call	printf
	
	pop	rbp		; restore stack
	mov	rax,0		; normal, no error, return value
	ret			; return

	; end main, code for function mat_vec
	
	section .bss
A:	resq 1			; address will be in rax
x:	resq 1			; address will be in rdx
y:	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]
n:	resq 1			; will be in r10
n8:	resq 1			; will be in r11 n*8
ij:	resq 1			; will be in r12 i*n+j
tmp:	resq 1
	
	section .data
zero:	dq 0.0
eight:	dq 8
dfmt1:	db "mat_vec n=%ld", 10, 0	; just to know execution got here
dfmtn8:	db "mat_vec n8=%ld", 10, 0	; debug
dfmti:	db "mat_vec i=%ld", 10, 0	; debug
dfmtj:	db "mat_vec j=%ld", 10, 0	; debug
dfmta:	db "A[i][j]=%f", 10, 0
dfmtx:	db "x[%1ld]=%f", 10, 0
dfmty:	db "y[%1ld]=%f", 10, 0
dfmtjf:	db "j loop finished", 10, 0
	
	section .text
	global mat_vec		; mat_vec(n,A,x,y)
mat_vec:
        push	rbp		; save rbp

	mov	[n],rdi		; value of n
	mov	[A],rsi		; save address of A matrix
	mov	[x],rdx		; save address of x vector
	mov	[y],rcx		; save address of y result
	
	mov	rdi, dfmt1
	mov	rsi, [n]	; print value of n
	mov	rax, 0
	call	printf

	mov	rax, [n]
	imul	qword [eight]
        mov	r11, rax	;  [*][]
	mov 	[n8], rax
	mov	r10, [n]	; <n
	
	mov	rdi, dfmtn8
	mov	rsi, [n8]	; print value of n*8
	mov	rax, 0
	call 	printf

	mov	rax, [A]	; address of A matrix
	mov	rdx, [x]	; address of x vector
	mov	rcx, [y]	; address of y vector
	
	mov	r8, 0		; i = 0
iloop:				; for(i=0; i<n; i++)
	mov	[i],r8
	mov	rdi, dfmti
	mov	rsi, [i]	; print value of i
	mov	rax, 0
	call	printf

	mov	rax, [A]	; address of A matrix
	mov	rdx, [x]	; address of x vector
	mov	rcx, [y]	; address of y vector
	mov	r11, [n8]	; n*8
	mov	r10, [n]	; n
	mov	r9, [j]		; j
	mov	r8, [i]		; i
	
	fld	qword [zero]
	fstp	qword [rcx+r8]  ; y[i] = 0.0

	mov	r9, 0		; j = 0
jloop:				; for(j=0; j<n; j++)
				; y[i] = y[i] + A[j][i]*x[j];
	mov	[j], r9
	mov	rdi, dfmtj
	mov	rsi, [j]	; print value of j
	mov	rax, 0
	call	printf

	mov	r9,[j]
	mov	rdx, [x]	; address of x matrix
	fld	qword [rdx+r9]  ; x[j]
	fstp	qword [tmp]
	mov	rdi, dfmtx
	mov	rsi, [j]
	movq	xmm0, qword [tmp]
	mov	rax, 1
	call	printf

	mov	r8, [i]
	mov	rcx, [y]	; address of y matrix
	fld	qword [rcx+r8]   ; y[i]
	fstp	qword [tmp]
	mov	rdi, dfmty
	mov	rsi, [i]
	movq	xmm0, qword [tmp]
	mov	rax, 1
	call	printf

	mov	rax, [i]	; i
	imul	qword [n]	; n
	add	rax, [j]	; j
	mov	r12,rax		; i*n+j  i*8*n+j*8 bytes
        mov	[ij],rax
	
	mov	rax, [A]	; address of A[i][j] 
	fld	qword [rax+r12]
	fstp	qword [tmp]
	mov	rdi, dfmta
	movq	xmm0, qword [tmp]
	mov	rax, 1
	call	printf
	
	mov	rax, [A]	; address of A matrix
	mov	rdx, [x]	; address of x vector
	mov	rcx, [y]	; address of y vector
	mov	r12, [ij]	; i*n+j
	mov	r11, [n8]	; n*8
	mov	r10, [n]	; n
	mov	r9, [j]		; j
	mov	r8, [i]		; i

	fld     qword [rax+r12]     ; A[j][i]
	fmul	qword [rdx+r9]	    ; x[j]
	fadd	qword [rcx+r8] 	    ; y[i]
	fstp	qword [rcx+r8]	    ; y[i]

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

	mov	rdi, dfmtjf
	mov	rax, 0
	call 	printf

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