x86-64 TUTORIAL: HILBERT MATRIX


LABOUCHERE SYSTEM PROGRAM USING x86-64 REGISTERS  |  NEW LIFE WITH JEKYLL!

The aim of solving this problem is to learn how to use the XMM registers for multiplication of floating point numbers. Matrix multiplication is a slow calculation especially if the floating point unit is used, and hence doing packed floating point calculations (if double precision is not required) might just be much faster. So this program will test that.

The Hilbert matrix is a simple matrix whose each element B(i, j) = 1/(i+j-1). For more details refer the Wikipedia Link.

The Hilbert matrix is a square matrix, so first we will take the size of the matrix as input from the user and then allocate the memory block for the matrix and store the elements of the matrix. Then, we shall calculate the square of the matrix and print it on screen. The floating point precision level that will be used is single precision or 4 bytes per number.

%macro prologue 0
        push    rbp
        mov     rbp,rsp
        push    rbx
        push    r12
        push    r13
        push    r14
        push    r15
%endmacro
%macro epilogue 0
        pop     r15
        pop     r14
        pop     r13
        pop     r12
        pop     rbx
        leave
        ret
%endmacro
;%define DEBUG_PRINT 1

section .rodata

    int_prompt      db "%d",0
    float_prompt    db "%0.03f | ",0
    prompt2         db "The square of the matrix is: ",10,0
    prompt1         db "Enter the size of the Hilbert matrix:",0

section .text
    global  main, create_hilbert_matrix, square_float_matrix, print_float_matrix
    extern  printf, scanf, malloc, free
    extern  print_int,print_nl
    main:
        push    rbp
        mov     rbp, rsp
        sub     rsp, 8
        push    rbx
        push    r12
        push    r13
        push    r14
        push    r15
        ; Enter the size of the matrix
        xor     rax, rax
        mov     rdi, dword prompt1
        call    printf
        xor     rax, rax
        lea     rsi, [rbp-8]
        mov     rdi, dword int_prompt
        call    scanf
        
        ; print the number entered    
        xor     rax, rax
        xor     rdi, rdi
        mov     rdi, [rbp-8]
        call print_int
        call print_nl
    
        ; Create a Hilbert Matrix
        xor     rax, rax
        xor     rdi, rdi
        mov     rdi, [rbp-8]
        call    create_hilbert_matrix
        mov     rbx, rax              ; store a pointer to the Hilbert matrix      
        push rbx
        push rax
        mov     rsi, rax
        xor     rdi, rdi
        mov     rdi, [rbp-8]
        call    print_float_matrix
        call print_nl
        pop rax 
        
        ; Calculate the Square of a floating point matrix
        ; Arguments are the size of the matrix and the pointer to the input matrix
        ; Return value is the pointer to the square of the input matrix
        mov     rsi, rax
        xor     rdi, rdi
        mov     rdi, [rbp-8]
        call    square_float_matrix
        mov     r12, rax              ; store a pointer to the square of the Hilbert matrix
        push r12
        ; Print the square of the Hilbert Matrix on the screen.
        mov     rsi, rax
        xor     rdi, rdi
        mov     rdi, [rbp-8]
        call    print_float_matrix
        call    print_nl
        pop r12
        ; End the program and free memory
        mov     rdi, r12
        xor     rax, rax
        call    free
        pop rbx
        mov     rdi, rbx
        xor     rax, rax
        call    free
        epilogue

    create_hilbert_matrix:
        prologue
        ; size of the matrix is in RDI
        mov     r15, rdi
        imul    rdi, rdi
        imul    rdi, 4          ; sizeof(float)
        ; size of the memory to be allocated in RDI now
        xor     rax, rax
        call    malloc
        ; store a copy of the memory location pointing to the matrix.
        mov     r11, rax
        mov     rbx, rax
        push rax
        push r11
        xor     rcx, rcx
outer_loop_start:
        cmp     rcx, r15
        jge     outer_loop_end
        xor     r12, r12
inner_loop_start:
        cmp     r12, r15
        jge     inner_loop_end

        ; each element is RCX+RDX+1
        mov     rdi, 1
        add     rdi, rcx
        adc     rdi, r12
      
        ; convert the integer in RDI to a single precision floating point in XMM0
        pxor xmm0, xmm0
        pxor xmm1, xmm1
        cvtsi2ss xmm0, rdi
        ; calculate the reciprocal of the single precision floating point in XMM0 and place it in XMM1
        rcpss   xmm1, xmm0
        ; place the value in XMM1 in memory allocated for the matrix
        movss    dword[rbx], xmm1

        ; Print the the Hilbert Matrix on the screen.
%ifdef DEBUG_PRINT
        push rbx
        push rax
        push rdx
        push rdi
        push rcx
        push r12
        push r15
        ; for printf we need to store the single precision as double precision
        cvtss2sd xmm0, dword[rbx]
        mov     rdi, dword float_prompt
        mov     rax, 1 ; Number of XMM registers to use
        call    printf
        pop r15
        pop r12
        pop rcx
        pop rdi
        pop rdx
        pop rax
        pop rbx
%endif

        add     rbx, 4          ; sizeof(float) 
        inc     r12
        jmp     inner_loop_start
inner_loop_end:
        ; Print the the Hilbert Matrix on the screen.
%ifdef DEBUG_PRINT
        push rbx
        push rax
        push rdx
        push rdi
        push rcx
        push r12
        push r15
        call print_nl
        pop r15
        pop r12
        pop rcx
        pop rdi
        pop rdx
        pop rax
        pop rbx
%endif
        inc     rcx
        jmp     outer_loop_start
outer_loop_end:
        pop     r11
        pop     rax
        mov     rax, r11
        epilogue

square_float_matrix:
         prologue
        ; the Hilbert matrix pointer is in RSI
        mov     r14, rsi
        ; size of the matrix is in RDI
        mov     r15, rdi
        imul    rdi, rdi
        imul    rdi, 4          ; sizeof(float)
        push   r15
        push   r14
        ; size of the memory to be allocated in RDI now
        xor     rax, rax
        call    malloc
        ; store a copy of the memory location pointing to the matrix.
        mov     rbx, rax
        mov     r12, rax
        pop     r14
        pop     r15
        push    rax
        push    rbx
        ; start squaring the matrix
        xor     rcx, rcx
forloop_1:
        cmp     rcx, r15
        jge     end_forloop_1
        xor     rsi, rsi
forloop_2:
        cmp     rsi, r15
        jge     end_forloop_2
        xor     rdi, rdi
        ; B[i][j] = 0.0 in XMM0
        pxor     xmm0, xmm0
        cvtsi2ss xmm0, rdi
forloop_3:
        cmp     rdi, r15
        jge     end_forloop_3
        ; B[i][j]+= A[i][k]*A[k][j]
        mov     rax, r15
        imul    rax, rcx
        add     rax, rdi
        mov     r8, r15
        imul    r8, rdi
        add     r8, rsi
        pxor    xmm1, xmm1
        pxor    xmm2, xmm2
        movss    xmm1, dword[r14 + rax*4]
        movss    xmm2, dword[r14 + r8*4]
        mulss   xmm1, xmm2
        addss   xmm0, xmm1
        inc     rdi
        jmp     forloop_3
end_forloop_3:
        ; B[i][j] = XMM0
        movss    dword[r12], xmm0
        add     r12, 4 ; sizeof(float)
        inc     rsi
        jmp     forloop_2
end_forloop_2:
        inc     rcx
        jmp     forloop_1
end_forloop_1:
        pop     rax
        pop     rbx
        mov     rax, rbx
        epilogue


print_float_matrix:
        prologue
        ; size of the matrix is in RDI and pointer is in RSI
        mov     r15, rdi
        mov     rbx, rsi
        xor     rcx, rcx
floop_1:
        cmp     rcx, r15
        jge      end_floop_1 
        xor     r12, r12
floop_2:
        cmp     r12, r15
        jge     end_floop_2

        push rbx
        push rax
        push rdx
        push rdi
        push rcx
        push r12
        push r15
        ; for printf we need to store the single precision as double precision
        cvtss2sd xmm0, dword[rbx]
        mov     rdi, dword float_prompt
        mov     rax, 1 ; Number of XMM registers to use
        call    printf
        pop r15
        pop r12
        pop rcx
        pop rdi
        pop rdx
        pop rax
        pop rbx

        add     rbx, 4
        inc     r12
        jmp     floop_2
end_floop_2:
        ; Print the the Hilbert Matrix on the screen.
        push rbx
        push rax
        push rdx
        push rdi
        push rcx
        push r12
        push r15
        call print_nl
        pop r15
        pop r12
        pop rcx
        pop rdi
        pop rdx
        pop rax
        pop rbx

        inc     rcx
        jmp     floop_1
end_floop_1:
        xor     rax, rax
        epilogue

If the above assembly program was called hilbert.asm then executing the following statements would create the executable hilbert.out.

$ yasm -f elf64 hilbert.asm
$ ld -m elf_x86_64 -dynamic-linker /lib64/ld-linux-x86-64.so.2 \
      /usr/lib/x86_64-linux-gnu/crt1.o /usr/lib/x86_64-linux-gnu/crti.o \
      hilbert.o asm_io.o /usr/lib/x86_64-linux-gnu/crtn.o -lc -o hilbert.out

The output is as follows for a 4x4 Hilbert matrix.

$ ./hilbert.out 
Enter the size of the Hilbert matrix:4
4
1.000 | 0.500 | 0.333 | 0.250 | 
0.500 | 0.333 | 0.250 | 0.200 | 
0.333 | 0.250 | 0.200 | 0.167 | 
0.250 | 0.200 | 0.167 | 0.143 | 

1.423 | 0.800 | 0.566 | 0.441 | 
0.800 | 0.463 | 0.333 | 0.262 | 
0.566 | 0.333 | 0.241 | 0.190 | 
0.441 | 0.262 | 0.190 | 0.151 | 

Download hilbert.asm, asm_io.inc and asm_io.asm.


LABOUCHERE SYSTEM PROGRAM USING x86-64 REGISTERS  |  NEW LIFE WITH JEKYLL!
SUPPORT THIS SITE
Donate DOGECOIN to DBevjMg3fd8C5oxZbV8sFpAffo6Tas1s8Q. DBevjMg3fd8C5oxZbV8sFpAffo6Tas1s8Q Donate BITCOIN to 19hrWWw1dPvBE1wVPfCnH8LqnUwsT3NsHW. 19hrWWw1dPvBE1wVPfCnH8LqnUwsT3NsHW
As an Amazon Associate I earn from qualifying purchases.