Я думаю, что лучший способ пойти, чтобы проверить как 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() "O(n^3)")
гораздо выше, чем одна из замен sgetrs()/dgetrs() "O(n^2.n_{rhs})")
. Поскольку максимальное число итераций в наборе 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 итераций