2016-03-26 7 views
1

В документации LAPACK говорится, что DSGESV (или ZCGESV для комплексных чисел) является:Когда использовать dsgesv против dgesv решить систему линейных уравнений

dsgesv и zcgesv смешивают точность итеративного уточнения подпрограмм для быстрого использования аппаратуры с одной точностью. Они сначала пытаются разложить матрицу в одиночной точности (dsgesv) или одиночную сложную точность (zcgesv) и использовать эту факторизацию в процедуре итеративной обработки для получения решения с двойной точностью 0 × 355 ×(dsgesv)/double complex precision (zcgesv) по норме качество обратной ошибки (см. ниже). Если подход завершается неудачно, метод переключается на факторизацию с двойной точностью или двойной комплексной точностью и вычисляет решение.

Итеративная утонченность не будет выигрышной стратегией, если коэффициент одиночной точности по сравнению с двойной точностью слишком мал. Разумная стратегия должна учитывать количество правых сторон и размер матрицы. Это может сделать с призывом к ilaenv в будущем. В настоящее время реализована итеративная обработка .

Но как я могу узнать, что такое отношение производительности одной точности к двойной точности? Существует предположение о том, чтобы учесть размер матрицы, но я не вижу, как именно размер матрицы приводит к оценке этого коэффициента производительности.

Может ли кто-нибудь прояснить эти вещи?

ответ

1

Я думаю, что лучший способ пойти, чтобы проверить как dgesv() и dsgesv() ...

Глядя на исходный код функции dsgesv() из Lapack, вот что dsgesv() пытается выполнить:

  • Cast матрицу A плавать As
  • вызов sgetrf(): LU факторизация, одинарной точности
  • Так LVE систему As.x=b с помощью LU факторизация по телефону sgetrs()
  • Compute двойной точности остаток r=b-Ax и решить As.x'=r с помощью sgetrs() снова, добавьте x=x+x'.

Последний шаг повторяется до тех пор, пока не будет достигнута двойная точность (максимум 30 итераций). Критерии определения успеха:

где является точность чисел с плавающей запятой двойной Precison (примерно 1e-13) и является размер матрицы. Если он терпит неудачу, dsgesv() возобновляет до dgesv(), так как он вызывает dgetrf() (факторизация), а затем dgetrs(). Следовательно, dsgesv() - алгоритм смешанной точности. См. Например, this article.

Наконец, dsgesv(), как ожидается, опережать dgesv() для небольшого числа правых частей и больших матриц, то есть, когда стоимость факторизации sgetrf()/dgetrf() гораздо выше, чем одна из замен sgetrs()/dgetrs(). Поскольку максимальное число итераций в наборе dsgesv() 30, приблизительный предел будет

Кроме того, sgetrf() должен proove значительно быстрее, чем dgetrf(). sgetrf() может быть быстрее из-за ограниченной доступной пропускной способности памяти или векторной обработки (смотрите SIMD, пример из SSE: инструкция ADDPS).

Аргумент iter из dsgesv() может быть протестирован для проверки полезности итеративной обработки. Если он отрицательный, итеративная обработка не удалась и использование dsgesv() было просто пустой тратой времени!

Здесь код С для сравнения и время dgesv(), sgesv(), dsgesv(). Он может быть скомпилирован gcc main.c -o main -llapacke -llapack -lblas Не стесняйтесь протестировать свою собственную матрицу!

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

#include <lapacke.h> 

int main(void){ 

    srand (time(NULL)); 

    //size of the matrix 
    int n=2000; 
    // number of right-hand size 
    int nb=3; 

    int nbrun=1000*100*100/n/n; 

    //memory initialization 
    double *aaa=malloc(n*n*sizeof(double)); 
    if(aaa==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    double *aa=malloc(n*n*sizeof(double)); 
    if(aa==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    double *bbb=malloc(n*nb*sizeof(double)); 
    if(bbb==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    double *x=malloc(n*nb*sizeof(double)); 
    if(x==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    double *bb=malloc(n*nb*sizeof(double)); 
    if(bb==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 

    float *aaas=malloc(n*n*sizeof(float)); 
    if(aaas==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    float *aas=malloc(n*n*sizeof(float)); 
    if(aas==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    float *bbbs=malloc(n*n*sizeof(float)); 
    if(bbbs==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    float *bbs=malloc(n*nb*sizeof(float)); 
    if(bbs==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 

    int *ipiv=malloc(n*nb*sizeof(int)); 
    if(ipiv==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 

    int i,j; 

    //matrix initialization 
    double cond=1e3; 
    for(i=0;i<n;i++){ 
     for(j=0;j<n;j++){ 
      if(j==i){ 
       aaa[i*n+j]=pow(cond,(i+1)/(double)n); 
      }else{ 
       aaa[i*n+j]=1.9*(rand()/(double)RAND_MAX-0.5)*pow(cond,(i+1)/(double)n)/(double)n; 
       //aaa[i*n+j]=(rand()/(double)RAND_MAX-0.5)/(double)n; 
       //aaa[i*n+j]=0; 
      } 
     } 
     bbb[i]=i; 
    } 

    for(i=0;i<n;i++){ 
     for(j=0;j<n;j++){ 
      aaas[i*n+j]=aaa[i*n+j]; 
     } 
     bbbs[i]=bbb[i]; 
    } 

    int k=0; 
    int ierr; 


    //estimating the condition number of the matrix 
    memcpy(aa,aaa,n*n*sizeof(double)); 
    double anorm; 
    double rcond; 
    //anorm=LAPACKE_dlange(LAPACK_ROW_MAJOR, 'i', n,n, aa, n); 
    double work[n]; 
    anorm=LAPACKE_dlange_work(LAPACK_ROW_MAJOR, 'i', n, n, aa, n, work); 
    ierr=LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n,aa, n,ipiv); 
    if(ierr<0){LAPACKE_xerbla("LAPACKE_dgetrf", ierr);} 
    ierr=LAPACKE_dgecon(LAPACK_ROW_MAJOR, 'i', n,aa, n,anorm,&rcond); 
    if(ierr<0){LAPACKE_xerbla("LAPACKE_dgecon", ierr);} 
    printf("condition number is %g\n",anorm,1./rcond); 

    //testing dgesv() 
    clock_t t; 
    t = clock(); 
    for(k=0;k<nbrun;k++){ 

     memcpy(bb,bbb,n*nb*sizeof(double)); 
     memcpy(aa,aaa,n*n*sizeof(double)); 



     ierr=LAPACKE_dgesv(LAPACK_ROW_MAJOR,n,nb,aa,n,ipiv,bb,nb); 
     if(ierr<0){LAPACKE_xerbla("LAPACKE_dgesv", ierr);} 

    } 

    //testing sgesv() 
    t = clock() - t; 
    printf ("dgesv()x%d took me %d clicks (%f seconds).\n",nbrun,t,((float)t)/CLOCKS_PER_SEC); 

    t = clock(); 
    for(k=0;k<nbrun;k++){ 

     memcpy(bbs,bbbs,n*nb*sizeof(float)); 
     memcpy(aas,aaas,n*n*sizeof(float)); 



     ierr=LAPACKE_sgesv(LAPACK_ROW_MAJOR,n,nb,aas,n,ipiv,bbs,nb); 
     if(ierr<0){LAPACKE_xerbla("LAPACKE_sgesv", ierr);} 

    } 

    //testing dsgesv() 
    t = clock() - t; 
    printf ("sgesv()x%d took me %d clicks (%f seconds).\n",nbrun,t,((float)t)/CLOCKS_PER_SEC); 

    int iter; 
    t = clock(); 
    for(k=0;k<nbrun;k++){ 

     memcpy(bb,bbb,n*nb*sizeof(double)); 
     memcpy(aa,aaa,n*n*sizeof(double)); 


     ierr=LAPACKE_dsgesv(LAPACK_ROW_MAJOR,n,nb,aa,n,ipiv,bb,nb,x,nb,&iter); 
     if(ierr<0){LAPACKE_xerbla("LAPACKE_dsgesv", ierr);} 

    } 
    t = clock() - t; 
    printf ("dsgesv()x%d took me %d clicks (%f seconds).\n",nbrun,t,((float)t)/CLOCKS_PER_SEC); 

    if(iter>0){ 
     printf("iterative refinement has succeded, %d iterations\n"); 
    }else{ 
     printf("iterative refinement has failed due to"); 
     if(iter==-1){ 
      printf(" implementation- or machine-specific reasons\n"); 
     } 
     if(iter==-2){ 
      printf(" overflow in iterations\n"); 
     } 
     if(iter==-3){ 
      printf(" failure of single precision factorization sgetrf() (ill-conditionned?)\n"); 
     } 
     if(iter==-31){ 
      printf(" max number of iterations\n"); 
     } 
    } 
    free(aaa); 
    free(aa); 
    free(bbb); 
    free(bb); 
    free(x); 


    free(aaas); 
    free(aas); 
    free(bbbs); 
    free(bbs); 

    free(ipiv); 

    return 0; 
} 

Выход для п = 2000:

число обусловленности 1475,26

dgesv() x2 взял меня 5260000 щелчки (5.260000 секунд).

sgesv() x2 взял меня 3560000 кликов (3,560000 секунд).

dsgesv() x2 взял меня за 3790000 кликов (3.790000 секунд).

итерационный уточнение преуспели, 11 итераций