VPF::Умножение двух матриц рекурсивным путём

ressac

Дата 2.11.2009, 19:16 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 345
Регистрация: 25.11.2006

Репутация: нет
Всего: 1

ну вот

, есть две матрицы и надо сделать рекурсивный алгоритм умножения, я что-то набросал, но запутался.

Anikmar

Дата 2.11.2009, 19:22 (ссылка)

(голосов:4)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2513
Регистрация: 26.11.2006
Где: Санкт-Петербург

Репутация: 9
Всего: 59

Цитата(ressac @ 2.11.2009, 19:16

)

я что-то набросал, но запутался.

Шпион однако. Невидимыми чернилами?

ressac

Дата 2.11.2009, 23:01 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 345
Регистрация: 25.11.2006

Репутация: нет
Всего: 1

да нет, я просто и не стал постить, просто чтоб стыдно не было.... ну вот если хотите, я его ещё не закончил

Код


#define N 10

void productoMatrix(int[N][N],int[N][N],int[N][N],int,int);
void productoFila(int[],int[N][N],int[],int,int);
void productoFilaAux(int[],int[N][N],int[],int,int);

int main()
{
    int m1[N][N]=
    {
        {1,2,3,4},
        {2,1,5,6},
        {3,5,1,7},
        {4,6,7,1}
    };

    int m2[N][N]=
    {
        {1,2,3,4},
        {2,1,5,6},
        {3,5,1,7},
        {4,6,7,1}
    };

    int mR[N][N]={0};

    int TAM = 4;
    int x,y;

    productoMatrix(m1,TAM-1,TAM-1
                   ,m2,TAM-1,TAM-1
                   ,mR,TAM-1,TAM-1);

    for (x=0;x<TAM;x++)
    {
        for (y=0;y<TAM;y++)
            printf("%4i",mR[x][y]);

        printf("\n");
    }
}

void productoMatrix(int m1[N][N],int fM1,int cM1,
                    int m2[N][N],int fM2,int cM2,
                    int mR[N][N],int fR,int cR)
{
    if (0 <= fM1)
    {
        productoMatrix(m1,fM1-1,cM1,
                       m2,fM2,cM2-1,
                       mR,fR-1,cR);

        productoFila(&m1[fMi],
                     m2,fM2,cM2,
                     &mR[fR],cR);
    }
}

void productoFila(int m1[],
                  int m2[N][N],int fM2,int cM2,
                  int mR[],int cR)
{
    if (0 <= cR)
    {
        productoFila(m1,
                     m2,fM2,cM2,
                     mR,cR-1);

        mR[cR] = productoFilaAux(m1,
                                 m2,cM2,
                                 fcComun);
    }
    else
        mR[cR] = 0;
}

int productoFilaAux(int m1[],
                    int m2[N][N],int cM2
                    int fcComun)
{
    int r;

    if (0 <= fcComun)
    {
        r += productoFilaAux(m1,
                             m2[N][N],
                             fM2,fcComun-1)

             +

             m1[fcComun] * m2[fcComun][cM2];
    }
    else
        r = 0;

    return r;
}

sdukshis

Дата 3.11.2009, 00:34 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 92
Регистрация: 23.3.2009

Репутация: нет
Всего: 1

Могу предложить черновой вариант алгоритма Штрассена.
Заранее прошу прощения за отсутствие комментариев

Код


#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;
}

ressac

Дата 3.11.2009, 00:54 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 345
Регистрация: 25.11.2006

Репутация: нет
Всего: 1

sdukshis, разве он рекурсивный?

Anikmar

Дата 3.11.2009, 01:08 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2513
Регистрация: 26.11.2006
Где: Санкт-Петербург

Репутация: 9
Всего: 59

Цитата(ressac @ 2.11.2009, 23:01

)

да нет, я просто и не стал постить, просто чтоб стыдно не было.... ну вот если хотите, я его ещё не закончил

Я наивно пердположил увидеть участок кода, в котором вы запутались smile

ressac

Дата 3.11.2009, 02:51 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 345
Регистрация: 25.11.2006

Репутация: нет
Всего: 1

сделал, с горем пополам smile

вроде оптимально, но я думаю можно ещё оптимальней,

может кто знает другой способ?

Код


#define N 10

void productoMatrix(int[N][N],int,
                    int[N][N],int,
                    int[N][N],int);

void productoFila(int[],
                  int[N][N],int,
                  int[],int);

int productoFilaAux(int[],
                    int [N][N],int,
                    int);

int main()
{
    int m1[N][N]=
    {
        {1,2,3,4},
        {2,1,5,6},
        {3,5,1,7},
        {4,6,7,1}
    };

    int m2[N][N]=
    {
        {1,2,3,4},
        {2,1,5,6},
        {3,5,1,7},
        {4,6,7,1}
    };

    int mR[N][N]={0};

    int TAM = 4;
    int x,y;

    productoMatrix(m1,TAM-1,
                   m2,TAM-1,
                   mR,TAM-1);

    for (x=0;x<TAM;x++)
    {
        for (y=0;y<TAM;y++)
            printf("%4i",mR[x][y]);

        printf("\n");
    }
}

void productoMatrix(int m1[N][N],int fM1,
                    int m2[N][N],int cM2,
                    int mR[N][N],int fcComun)
{
    if (0 <= fM1)
    {
        productoMatrix(m1,fM1-1,
                       m2,cM2,
                       mR,fcComun);

        productoFila(&m1[fM1],
                     m2,cM2,
                     &mR[fM1],fcComun);
    }
}

void productoFila(int m1[],
                  int m2[N][N],int cM2,
                  int mR[],int fcComun)
{
    if (0 <= cM2)
    {
        productoFila(m1,
                     m2,cM2-1,
                     mR,fcComun);

        mR[cM2] = productoFilaAux(m1,
                                  m2,cM2,
                                  fcComun);
    }
    else
        mR[cM2+1] = 0;
}

int productoFilaAux(int m1[],
                    int m2[N][N],int cM2,
                    int fcComun)
{
    int r;

    if (0 <= fcComun)
    {
        r = productoFilaAux(m1,
                            m2,cM2,
                            fcComun-1)

            +

            m1[fcComun] * m2[fcComun][cM2];
    }
    else
        r = 0;

    return r;
}

sdukshis

Дата 3.11.2009, 13:43 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 92
Регистрация: 23.3.2009

Репутация: нет
Всего: 1

Цитата(ressac @ 3.11.2009, 00:54)

sdukshis, разве он рекурсивный?

Конечно. Это хорошо видно в строках 157-163 shtras_mul() вызывает саму себя

FCM

Дата 3.11.2009, 20:57 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: нет
Всего: 9

Позвольте поинтересоваться, рекурсивное решение подобных задач - это академическое упражнение или насущная задача, имеющая особую практическую ценность в плане производительности?.

ressac

Дата 3.11.2009, 21:53 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 345
Регистрация: 25.11.2006

Репутация: нет
Всего: 1

упражнение

Anikmar

Дата 4.11.2009, 11:39 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Завсегдатай
Сообщений: 2513
Регистрация: 26.11.2006
Где: Санкт-Петербург

Репутация: 9
Всего: 59

Цитата(FCM @ 3.11.2009, 20:57

)

вольте поинтересоваться, рекурсивное решение подобных задач - это академическое упражнение или насущная задача, имеющая особую практическую ценность в плане производительности?.

А что, в общем виде рекурсивный метод решения дает ощутимый прирост производительности? Мне всегда казалось, что кроме простоты кода никаких преимуществ рекурсия не дает.

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| C/C++: Общие вопросы \| Следующая тема »