Преобразование симметричной матрицы между упакованным и полным хранилищем?

Question

Я новичок в числовой линейной алгебре, и я только начал использовать LAPACK и BLAS.

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

Я нашел dtrttp, который я могу использовать для преобразования полной симметричной матрицы с двойной точностью в упакованное хранилище. Однако эти подпрограммы предназначены для треугольных матриц, поэтому соответствующий dtpttr заполняет только треугольник полной матрицы. Как я могу заполнить вторую половину?

Answer 1

Очевидное решение состоит в том, чтобы симметризировать матрицу с помощью кода «home-made/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 эти петли занимают 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; 
}