在打包和完全存储之间转换对称矩阵？

KQS 发表于 Dev

157

品质管理系统

我是数值线性代数的新手，我刚刚开始使用LAPACK和BLAS。

是否有例程可以在打包存储和完整存储之间复制/转换对称矩阵？

我发现了dtrttp，可以用来将双精度全对称矩阵转换为打包存储。但是，这些例程适用于三角形矩阵，因此相应的函数dtpttr仅填充整个矩阵的三角形。我如何填补另一半？

弗朗西斯

显而易见的解决方案是通过“自制/ diy”代码对称化矩阵，这是重新发明轮子的风险。在之后编写对称矩阵所需的for循环非常容易dtpttr。

for(i=0;i<n;i++){
  for(j=i+1;j<n;j++){
    a[i*n+j]=a[j*n+i];
  }
}

它对您的应用程序足够有效吗？在10000x10000矩阵上，这些for循环在我的PC上dtpttr持续0.88s ，而持续0.24s。

这是测试代码。用以下命令编译gcc main.c -o main -llapack -lblas -lm：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

void dtrttp_(char* UPLO,int* N,double* A,int* LDA,double* AP,int* INFO);
void dtpttr_(char* UPLO,int* N,double* AP,double* A,int* LDA,int* INFO);
void daxpy_(int* N,double* DA,double* DX,int* INCX,double* DY,int* INCY);

void dtpttf_(char* TRANSR,char* UPLO,int* N,double* AP,double* ARF,int* INFO);

int main(int argc, char **argv)
{

    int n=10;
    int info;

    double *a=malloc(n*n*sizeof(double));
    double *packed=malloc((n*(n+1))/2*sizeof(double));

    int i,j;
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            a[i*n+j]=i+j;
        }
    }

    printf("matrix before pack\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf("%g ",a[i*n+j]);
        }
        printf("\n");
    }

    printf("\n");
    //pack
    dtrttp_("U",&n,a,&n,packed,&info);

    //unpack
    memset(a,0,n*n*sizeof(double));
    dtpttr_("U",&n,packed,a,&n,&info);

    for(i=0;i<n;i++){
        for(j=i+1;j<n;j++){
            a[i*n+j]=a[j*n+i];
        }
    }

    printf("matrix after unpack\n");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf("%g ",a[i*n+j]);
        }
        printf("\n");
    }

    free(a);
    free(packed);

    printf("timing...\n");

    n=10000;

    a=malloc(n*n*sizeof(double));
    packed=malloc((n*(n+1))/2*sizeof(double));

    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            a[i*n+j]=i+j;
        }
    }

    //pack
    dtrttp_("U",&n,a,&n,packed,&info);

    //unpack
    memset(a,0,n*n*sizeof(double));
    clock_t t;
    t = clock();
    dtpttr_("U",&n,packed,a,&n,&info);
    t = clock() - t;
    printf ("dtpttr took %f seconds.\n",((double)t)/CLOCKS_PER_SEC);
    t = clock();
    for(i=0;i<n;i++){
        for(j=i+1;j<n;j++){
            a[i*n+j]=a[j*n+i];
        }
    }
    t = clock() - t;
    printf ("symmetrize took %f seconds.\n",((double)t)/CLOCKS_PER_SEC);
    free(a);
    free(packed);

    return 0;
}

本文收集自互联网，转载请注明来源。

如有侵权，请联系 [email protected] 删除。