在LAPACK文档中，它声明DSGESV（或ZCGESV表示复数）是：
dsgesv和zcgesv是混合精度迭代求精
用于开发快速单精度硬件的子程序。他们先
尝试将矩阵分解为单精度（dsgesv）或单精度
复杂精度（zcgesv）并在
用迭代求精法求双解
标准精度（dsgesv）/双复精度（zcgesv）
向后误差质量（见下文）如果进近失败
切换到双精度或双复精度
分别进行因式分解和求解。
如果
单精度性能与双精度性能之比
太小了。一个合理的策略应该是
右手边和矩阵的大小考虑在内。这可能
以后再打个电话给ilanv目前，迭代
实现了优化。
但是，我怎么知道单精度性能与双精度性能的比率是多少呢？有人建议考虑矩阵的大小，但我不知道矩阵的大小是如何导致对这个性能比的估计的。
有人能澄清这些事情吗？
最佳答案
我想最好的办法是同时测试dgesv()和dsgesv()。
查看lapack函数dsgesv()的源代码，这里是dsgesv()试图执行的操作：
将矩阵A转换为浮点As
调用：lu因子分解，单精度
通过调用lu因子分解来解决系统sgetrf()
计算双精度残差As.x=b并再次使用sgetrs()求解r=b-Ax，添加As.x'=r。
重复最后一步，直到达到双精度（最多30次迭代）。定义成功的标准是：
双精度浮点数（大约1E-13）的精度在哪里，是矩阵的大小。如果失败，sgetrs()将恢复到x=x+x'，因为它调用dsgesv()（因子分解），然后dgesv()。因此dgetrf()是一种混合精度算法。例如，请参见this article。
最后，对于少量的右手边和大型矩阵，dgetrs()的性能有望优于dsgesv()，也就是说，当因子分解的成本远高于替代的成本时。由于dsgesv()中迭代集的最大数目为30，近似极限将是
此外，dgesv()必须明显快于sgetrf()。dgetrf()由于可用内存带宽或矢量处理有限，因此速度可能更快（请查找来自SSE的SIMD示例：指令sgetrs()）。
可以测试dgetrs()的参数dsgesv()来检查迭代求精是否有用。如果是负数，则迭代求精失败，使用sgetrf()只是浪费时间！
这里有一个C代码来比较和计时它可以由dgetrf()编译，请随意测试您自己的矩阵！
#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;
}
n=2000时的输出：
条件号是1475.26
我点击了5260000次（5.260000秒）。
sgesv（）x2点击了3560000次（3.560000秒）。
dsgesv（）x2花费了我3790000次点击（3.790000秒）。
迭代优化成功，11次迭代
关于algorithm - 何时使用dsgesv与dgesv求解线性方程组，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/36241631/