main.cpp
#include <stdlib.h> #include <stdio.h> extern "C" double* MatrixMul(const double* m1, int nr1, int nc1, const double* m2, int nr2, int nc2); void MatrixPrint(const double* m, int nr, int nc, const char* s) { printf("%s\n", s); if (m != NULL) { for (int i = 0; i < nr; i++) { for (int j = 0; j < nc; j++) { double m_val = m[i * nc + j]; printf("%8.1lf ", m_val); } printf("\n"); } } else printf("NULL pointer\n"); } double* MatrixMulCpp(const double* m1, int nr1, int nc1, const double* m2, int nr2, int nc2) { if ((nr1 < 0) || (nc1 < 0) || (nr2 < 0) || (nc2 < 0)) return NULL; if (nc1 != nr2) return NULL; double* m3 = (double*)malloc(nr1 * nc2 * sizeof(double)); for (int i = 0; i < nr1; i++) { for (int j = 0; j < nc2; j++) { double sum = 0; for (int k = 0; k < nc1; k++) { double m1_val = m1[i * nc1 + k]; double m2_val = m2[k * nc2 + j]; sum += m1_val * m2_val; } m3[i * nc2 + j] = sum; } } return m3; } void MatrixMul1(void) { const int nr1 = 3; const int nc1 = 2; const int nr2 = 2; const int nc2 = 3; double m1[nr1 * nc1] = { 6, 2, 4, 3, -5, -2 }; double m2[nr2 * nc2] = { -2, 3, 4, -3, 6, 7 }; double* m3_a = MatrixMulCpp(m1, nr1, nc1, m2, nr2, nc2); double* m3_b = MatrixMul(m1, nr1, nc1, m2, nr2, nc2); printf("\nResults for MatrixMul1()\n"); MatrixPrint(m1, nr1, nc1, "Matrix m1"); MatrixPrint(m2, nr2, nc2, "Matrix m2"); MatrixPrint(m3_a, nr1, nc2, "Matrix m3_a"); MatrixPrint(m3_b, nr1, nc2, "Matrix m3_b"); free(m3_a); free(m3_b); } void MatrixMul2(void) { const int nr1 = 2; const int nc1 = 3; const int nr2 = 3; const int nc2 = 4; double m1[nr1 * nc1] = { 5, -3, 2, -2, 5, 4 }; double m2[nr2 * nc2] = { 7, -4, 3, 3, 2, 6, -2, 5, 4, 9, 3, 5 }; double* m3_a = MatrixMulCpp(m1, nr1, nc1, m2, nr2, nc2); double* m3_b = MatrixMul(m1, nr1, nc1, m2, nr2, nc2); printf("\nResults for MatrixMul2()\n"); MatrixPrint(m1, nr1, nc1, "Matrix m1"); MatrixPrint(m2, nr2, nc2, "Matrix m2"); MatrixPrint(m3_a, nr1, nc2, "Matrix m3_a"); MatrixPrint(m3_b, nr1, nc2, "Matrix m3_b"); free(m3_a); free(m3_b); } int main(int argc, char* argv[]) { MatrixMul1(); MatrixMul2(); return 0; }
matrixmul.asm
; Name: matrixmul.asm ; ; Build: g++ -c main.cpp -o main.o ; nasm -f elf64 -o matrixmul.o matrixmul.asm ; g++ -o matrixmul matrixmul.o main.o ; ; Source: Modern x86 Assembly Language Programming p.546 ; extern "C" double* MatrixMul(const double* m1, int nr1, int nc1, const double* m2, int nr2, int nc2); ; ; Description: The following function computes the product of two ; matrices. global MatrixMul extern malloc section .text MatrixMul: ; GCC Windows ; Registers: rdi m1 rcx ; rsi nr1 rdx ; rdx nc1 r8 ; rcx m2 r9 ; r8 nr2 stack ; r9 nc2 stack ; Create stackframe push rbp mov rbp,rsp sub rsp,8 ;align stack push r12 ;save non-volatile register push r13 push r14 push r15 push rbx ; assume error xor rax,rax ;return NULL pointer when error ; Verify the matrix size values. movsxd rsi,esi ;rsi = nr1 test rsi,rsi jle .done ;jump if nr1 <= 0 movsxd rdx,edx ;rdx = nc1 test rdx,rdx jle .done ;jump if nc1 <= 0 movsxd r8,r8d ;r8 = nr2 test r8,r8 jle .done ;jump if nr2 <= 0 movsxd r9,r9d ;r9 = nc2 test r9,r9 jle .done ;jump if nc2 <= 0 cmp rdx,r8 jne .done ;jump if nc1 != nr2 ; Allocate storage ; save m1, m2, nc1 and nc2 on stack, these are altered by malloc() push rdi push rsi push rcx push rdx push r9 ;calculate number of bytes of memory to allocate mov rdi,rsi ;rdi = nr1 imul rdi,r9 ;rdi = nr1 * nc2 shl rdi,3 ;rdi = nr1 * nc2 * size real8 call malloc mov rbx,rax ;rbx = ptr to m3 ;restore m1, m2, nc1 and nc2 from stack pop r9 pop rdx pop rcx pop rsi pop rdi ; Initialize row index i xor r14,r14 ;i = 0 ; Initialize column index j .lp1: xor r15,r15 ;j = 0 ; Initialize sum and index k .lp2: xorpd xmm4,xmm4 ;sum = 0; xor r10,r10 ;k = 0; ; Calculate sum += m1[i * nc1 + k] * m2[k * nc2 + j] .lp3: mov rax,r14 ;rax = i imul rax,rdx ;rax = i * nc1 add rax,r10 ;rax = i * nc1 + k movsd xmm0,[rdi+rax*8] ;xmm0 = m1[i * nc1 + k] mov r11,r10 ;r11 = k; imul r11,r9 ;r11 = k * nc2 add r11,r15 ;r11 = k * nc2 + j movsd xmm1,[rcx+r11*8] ;xmm1 = m2[k * nc2 + j] mulsd xmm0,xmm1 ;xmm0 = m1[i * nc1 + k] * m2[k * nc2 + j] addsd xmm4,xmm0 ;update sum inc r10 ;k++ cmp r10,rdx jl .lp3 ;jump if k < nc1 ; Save sum to m3[i * nc2 + j] mov rax,r14 ;rax = i imul rax,r9 ;rax = i * nc2 add rax,r15 ;rax = i * nc2 + j movsd qword[rbx+rax*8],xmm4 ;m3[i * nc2 + j] = sum ; Update loop counters and repeat until done inc r15 ;j++ cmp r9,r15 jl .lp2 ;jump if j < nc2 inc r14 ;i++ cmp r14,rsi jl .lp1 ;jump if i < nr1 mov rax,rbx ;rax = ptr to m3 .done: pop rbx ;restore non-volatile registers pop r15 pop r14 pop r13 pop r12 mov rsp,rbp ;restore stack pop rbp ret
build
g++ -c main.cpp -o main.o nasm -f elf64 -o matrixmul.o matrixmul.asm g++ -o matrixmul matrixmul.o main.o