#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mkl.h>
#include "Prototipos.h"
/* Resolver en SECUENCIAL C = A*B
Entrada:
A es una matriz de dimensiones N x K, almacenada por filas (FiL_A=N, Col_A=K, lead dimension lda=K)
B es una matriz de dimensiones K x M, almacenada por filas (FiL_B=K, Col_B=M, lead dimension ldb=M)
BT es la matriz transpuesta de B de dimensiones M x K, almacenada por filas (Col_B =M, Fil_B=K, lead dimension ldbt=K)
Salida
C es una matriz de dimensiones N x M, almacenada por filas (FiL_C=N, Col_C=M, lead dimension ldc=M)
*/
void sec_dgemm(int m, int n, int k, double alfa, double *A, int lda, double *B, int ldb, double beta, double *C, int ldc);
double* Transpuesta(int m, int n, double *M);
int main(int argc, char *argv[])
{
int
Fil_A, Fil_B, Fil_C,
Col_A, Col_B, Col_C,
seed;
double
*MatA = NULL,
*MatB = NULL,
*MatC = NULL,
*MatLIB = NULL,
*MatTB = NULL,
alfa = 1.0,
beta = 0.0,
elapsed, ucpu, scpu;
if (argc != 5) {
printf("Uso: %s <Filas C> <Columnas C> <Filas B> <semilla>\n", argv[0]);
return -1;
}
Fil_C = Fil_A = atoi(argv[1]);
Col_C = Col_B = atoi(argv[2]);
Fil_B = Col_A = atoi(argv[3]);
seed = atoi(argv[4]);
/* Comprobando que ninguna dimension es nula */
if ((Fil_A == 0) || (Col_A == 0) || (Col_B == 0)) {
printf("Error 1: Alguna dimension es nula.\n");
return -1;
}
/* Reservando adecuadamente memoria para las matrices */
/* Recuerde que son 4, MatA, MatB, MatC y MatLIB. Las */
/* dimensiones de MatLIB son las mismas que para MatC */
MatA = (double *)malloc((Fil_A*Col_A) * sizeof(double));
MatB = (double *)malloc((Fil_B*Col_B) * sizeof(double));
MatC = (double *)malloc((Fil_A*Col_B) * sizeof(double));
MatLIB = (double *)malloc((Fil_A*Col_B) * sizeof(double));
MatTB = (double *)malloc((Col_B*Fil_B) * sizeof(double));
/* Generando aleatoriamente las matrices MatA y MatB */
Genera_Mat(MatA, Fil_A, Col_A, seed);
Genera_Mat(MatB, Fil_B, Col_B, seed + (seed/3));
/* Resolviendo el problema con la MKL */
Ctimer(&elapsed, &ucpu, &scpu, 0);
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, Fil_A, Col_B, Col_A, alfa, MatA, Col_A, MatB, Col_B, beta, MatLIB, Col_C);
Ctimer(&elapsed, &ucpu, &scpu, 1);
printf("El tiempo empleado por MKL es %2.7E segundos.\n", elapsed);
/* Resolviendo el problema Matriz x Matriz usando Naive modificado */
Ctimer(&elapsed, &ucpu, &scpu, 0);
/* Llame aqui a la funcion disenada en la practica02. MatC sera la matriz resultante */
MatTB = Transpuesta(Fil_B, Col_B, MatB);
sec_dgemm(Fil_A, Col_B, Col_A, alfa, MatA, Col_A, MatTB, Fil_B, beta, MatC, Col_C);
Ctimer(&elapsed, &ucpu, &scpu, 1);
printf("El tiempo empleado por Naive es %2.7E segundos.\n", elapsed);
printf("\nError Forward MKL vs Normal: %8.5E\n", Norma1(Fil_C, Col_C, MatLIB, MatC));
free(MatA); MatA = NULL;
free(MatB); MatB = NULL;
free(MatC); MatC = NULL;
free(MatLIB); MatLIB = NULL;
free(MatTB); MatTB = NULL;
return 0;
}
void sec_dgemm(int m, int n, int k, double alfa, double *A, int lda, double *B, int ldb, double beta, double *C, int ldc){
/*Calculo matricial*/
int i, j ,r;
for(i=0; i<m; i++){
for(j=0; j<n; j++){
C[i*ldc+j] =0;
for(r=0; r<k; r++){
C[i*ldc+j] = C[i*ldc+j] + A[i*lda+r] * B[r*ldb+j];
}
}
}
}
double* Transpuesta(int m, int n, double *M){
/*Funcion que hace la transpuesta de una matriz*/
double *MAux= NULL;
MAux = (double *)malloc((n*m) * sizeof(double));
int i,j;
for(i=0; i<n; i++){
for(j=0; j<m; j++)
MAux[j*m+i] = M[i*m+j];
}
free(MAux);
return MAux;
}