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

Question

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

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

我找到了dtrttp，我可以使用它来将双精度全对称矩阵转换为压缩存储。但是，这些例程适用于三角矩阵，因此对应的dtpttr 仅填充整个矩阵的三角形。另一半怎么填？

Answer 1

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

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

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

这里是测试代码。用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;
}