Могу предложить черновой вариант алгоритма Штрассена.
Заранее прошу прощения за отсутствие комментариев
| Код |
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <string.h>
int matrix_mul(double *a,double *b, double *c,const size_t n);
int matrix_add(double *a,double *b, double *c,const size_t n);
int matrix_sub(double *a,double *b, double *c,const size_t n);
int shtras_mul(double *a,double *b, double *c,const size_t n);
int shtras_muln(double *a,double *b, double *c,const size_t n);
void matrix_print(double *a,const size_t n);
int main(){
const size_t n=1024;
double *a,*b,*c;
a=(double*)malloc(n*n*sizeof(double));
b=(double*)malloc(n*n*sizeof(double));
c=(double*)malloc(n*n*sizeof(double));
size_t i,j;
srand(time(NULL));
for(i=0;i!=n*n;++i){
a[i]=(double)(rand()%10);
b[i]=(double)(rand()%10);
}
clock_t starttime=clock();
shtras_mul(a,b,c,n);
printf("shtras_mul takes %d seconds\n",(clock()-starttime)/CLOCKS_PER_SEC);
starttime=clock();
matrix_mul(a,b,c,n);
printf("matrix_mul takes %d seconds\n",(clock()-starttime)/CLOCKS_PER_SEC);
/*
matrix_print(a,n);
printf("\n");
matrix_print(b,n);
printf("\n");
matrix_mul(a,b,c,n);
matrix_print(c,n);
printf("\n");
shtras_muln(a,b,c,n);
matrix_print(c,n);
printf("\n");
*/
// shtras_mul(a,b,c,n);
// shtras_muln(a,b,c,n);
free(a);
free(b);
free(c);
return EXIT_SUCCESS;
}
void matrix_print(double *a,const size_t n){
size_t i,j;
for(i=0;i!=n;++i){
for (j=0;j!=n;++j){
printf("%5g ",*(a+i*n+j));
}
printf("\n");
}
}
int matrix_mul(double *a,double *b, double *c,const size_t n){
size_t i,j,k;
double val;
for(i=0;i!=n;++i){
for(j=0; j!=n; ++j){
c[i*n+j] = 0.0;
}
for(k=0;k!=n;++k){
val=a[i*n+k];
for(j=0;j!=n;++j){
c[i*n+j]+=val*b[k*n+j];
}
}
}
return 0;
};
int matrix_add(double *a,double *b, double *c,const size_t n){
size_t i;
for(i=0;i!=n*n;++i)
c[i]=a[i]+b[i];
return 0;
}
int matrix_sub(double *a,double *b, double *c,const size_t n){
size_t i;
for(i=0;i!=n*n;++i)
c[i]=a[i]-b[i];
return 0;
}
int shtras_mul(double *a,double *b, double *c,const size_t n){
if(n <= 64){
matrix_mul(a,b,c,n);
return 0;
}
double *A1,*A2,*A3,*A4,*A5,*A6,*A7;
double *B1,*B2,*B3,*B4,*B5,*B6,*B7;
double *P1,*P2,*P3,*P4,*P5,*P6,*P7;
const size_t m = n/2, m2=m*m;
A1=(double*)malloc(m2*sizeof(double));
A2=(double*)malloc(m2*sizeof(double));
A3=(double*)malloc(m2*sizeof(double));
A4=(double*)malloc(m2*sizeof(double));
A5=(double*)malloc(m2*sizeof(double));
A6=(double*)malloc(m2*sizeof(double));
A7=(double*)malloc(m2*sizeof(double));
B1=(double*)malloc(m2*sizeof(double));
B2=(double*)malloc(m2*sizeof(double));
B3=(double*)malloc(m2*sizeof(double));
B4=(double*)malloc(m2*sizeof(double));
B5=(double*)malloc(m2*sizeof(double));
B6=(double*)malloc(m2*sizeof(double));
B7=(double*)malloc(m2*sizeof(double));
P1=(double*)malloc(m2*sizeof(double));
P2=(double*)malloc(m2*sizeof(double));
P3=(double*)malloc(m2*sizeof(double));
P4=(double*)malloc(m2*sizeof(double));
P5=(double*)malloc(m2*sizeof(double));
P6=(double*)malloc(m2*sizeof(double));
P7=(double*)malloc(m2*sizeof(double));
size_t i,j,k;
for(i=0;i!=m;++i){
for(j=0;j!=m;++j){
k=i*m+j;
A1[k]=a[i*n+j];
B1[k]=b[i*n+j+m]-b[(i+m)*n+j+m];
A2[k]=a[i*n+j]+a[i*n+j+m];
B2[k]=b[(i+m)*n+j+m];
A3[k]=a[(i+m)*n+j]+a[(i+m)*n+j+m];
B3[k]=b[i*n+j];
A4[k]=a[(i+m)*n+j+m];
B4[k]=b[(i+m)*n+j]-b[i*n+j];
A5[k]=a[i*n+j]+a[(i+m)*n+j+m];
B5[k]=b[i*n+j]+b[(i+m)*n+j+m];
A6[k]=a[i*n+j+m]-a[(i+m)*n+j+m];
B6[k]=b[(i+m)*n+j]+b[(i+m)*n+j+m];
A7[k]=a[i*n+j]-a[(i+m)*n+j];
B7[k]=b[i*n+j]+b[i*n+j+m];
}
}
shtras_mul(A1,B1,P1,m);
shtras_mul(A2,B2,P2,m);
shtras_mul(A3,B3,P3,m);
shtras_mul(A4,B4,P4,m);
shtras_mul(A5,B5,P5,m);
shtras_mul(A6,B6,P6,m);
shtras_mul(A7,B7,P7,m);
for(i=0;i!=m;++i){
for(j=0;j!=m;++j){
k=i*m+j;
c[i*n+j]=P5[k]+P4[k]-P2[k]+P6[k];
c[i*n+j+m]=P1[k]+P2[k];
c[(i+m)*n+j]=P3[k]+P4[k];
c[(i+m)*n+j+m]=P5[k]+P1[k]-P3[k]-P7[k];
}
}
free(A1);free(A2);free(A3);free(A4);free(A5);free(A6);free(A7);
free(B1);free(B2);free(B3);free(B4);free(B5);free(B6);free(B7);
free(P1);free(P2);free(P3);free(P4);free(P5);free(P6);free(P7);
return 0;
};
int shtras_muln(double *a,double *b, double *c,const size_t n){
size_t m = 2;
while(n>m) m<<=1;
if (n==m){
shtras_mul(a,b,c,n);
return 0;
}
double *A,*B,*C;
A=(double*)malloc(m*m*sizeof(double));
B=(double*)malloc(m*m*sizeof(double));
C=(double*)malloc(m*m*sizeof(double));
memset(A,0,m*m*sizeof(double));memset(B,0,m*m*sizeof(double));
size_t i,j;
for(i=0;i!=n;++i){
for(j=0;j!=n;++j){
A[i*m+j]=a[i*n+j];
B[i*m+j]=b[i*n+j];
}
}
for(i=n;i!=m;++i){
A[i*m+i]=1;
B[i*m+i]=1;
}
shtras_mul(A,B,C,m);
for(i=0;i!=n;++i){
for(j=0;j!=n;++j){
c[i*n+j]=C[i*m+j];
}
}
free(A);free(B);free(C);
return 0;
}
|
|
, есть две матрицы и надо сделать рекурсивный алгоритм умножения, я что-то набросал, но запутался. 
