在封裝和完整存儲之間轉換對稱矩陣？

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