linpackdoc文档格式.docx
《linpackdoc文档格式.docx》由会员分享,可在线阅读,更多相关《linpackdoc文档格式.docx(25页珍藏版)》请在冰豆网上搜索。
staticREALaa[200][200],a[200][201],b[200],x[200];
REALcray,ops,total,norma,normx;
REALresid,residn,eps,t1,tm,tm2;
REALepslon(),second(),kf;
staticintipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
lda=201;
ldaa=200;
cray=.056;
n=100;
fprintf(stdout,ROLLING);
fprintf(stdout,PREC);
fprintf(stdout,"
PrecisionLinpack\n\n"
);
fprintf(stderr,ROLLING);
fprintf(stderr,PREC);
fprintf(stderr,"
ops=(2.0e0*(n*n*n))/3.0+2.0*(n*n);
matgen(a,lda,n,b,&
norma);
t1=second();
dgefa(a,lda,n,ipvt,&
info);
time[0][0]=second()-t1;
dgesl(a,lda,n,ipvt,b,0);
time[1][0]=second()-t1;
total=time[0][0]+time[1][0];
/*computearesidualtoverifyresults.*/
for(i=0;
i<
n;
i++){
x[i]=b[i];
}
b[i]=-b[i];
dmxpy(n,b,n,lda,x,a);
resid=0.0;
normx=0.0;
resid=(resid>
fabs((double)b[i]))
?
resid:
fabs((double)b[i]);
normx=(normx>
fabs((double)x[i]))
normx:
fabs((double)x[i]);
eps=epslon((REAL)ONE);
residn=resid/(n*norma*normx*eps);
printf("
norm.residresidmachep"
printf("
x[0]-1x[n-1]-1\n"
%8.1f%16.8e%16.8e%16.8e%16.8e\n"
(double)residn,(double)resid,(double)eps,
(double)x[0]-1,(double)x[n-1]-1);
fprintf(stderr,"
timesarereportedformatricesoforder%5d\n"
n);
fprintf(stderr,"
dgefadgesltotalkflopsunit"
ratio\n"
time[2][0]=total;
time[3][0]=ops/(1.0e3*total);
time[4][0]=2.0e3/time[3][0];
time[5][0]=total/cray;
timesforarraywithleadingdimensionof%5d\n"
lda);
print_time(0);
time[0][1]=second()-t1;
time[1][1]=second()-t1;
total=time[0][1]+time[1][1];
time[2][1]=total;
time[3][1]=ops/(1.0e3*total);
time[4][1]=2.0e3/time[3][1];
time[5][1]=total/cray;
time[0][2]=second()-t1;
time[1][2]=second()-t1;
total=time[0][2]+time[1][2];
time[2][2]=total;
time[3][2]=ops/(1.0e3*total);
time[4][2]=2.0e3/time[3][2];
time[5][2]=total/cray;
ntimes=NTIMES;
tm2=0.0;
ntimes;
tm=second();
matgen(a,lda,n,b,&
tm2=tm2+second()-tm;
dgefa(a,lda,n,ipvt,&
time[0][3]=(second()-t1-tm2)/ntimes;
dgesl(a,lda,n,ipvt,b,0);
time[1][3]=(second()-t1)/ntimes;
total=time[0][3]+time[1][3];
time[2][3]=total;
time[3][3]=ops/(1.0e3*total);
time[4][3]=2.0e3/time[3][3];
time[5][3]=total/cray;
print_time
(1);
print_time
(2);
print_time(3);
matgen(aa,ldaa,n,b,&
dgefa(aa,ldaa,n,ipvt,&
time[0][4]=second()-t1;
dgesl(aa,ldaa,n,ipvt,b,0);
time[1][4]=second()-t1;
total=time[0][4]+time[1][4];
time[2][4]=total;
time[3][4]=ops/(1.0e3*total);
time[4][4]=2.0e3/time[3][4];
time[5][4]=total/cray;
time[0][5]=second()-t1;
time[1][5]=second()-t1;
total=time[0][5]+time[1][5];
time[2][5]=total;
time[3][5]=ops/(1.0e3*total);
time[4][5]=2.0e3/time[3][5];
time[5][5]=total/cray;
time[0][6]=second()-t1;
time[1][6]=second()-t1;
total=time[0][6]+time[1][6];
time[2][6]=total;
time[3][6]=ops/(1.0e3*total);
time[4][6]=2.0e3/time[3][6];
time[5][6]=total/cray;
tm2=0;
matgen(aa,ldaa,n,b,&
dgefa(aa,ldaa,n,ipvt,&
time[0][7]=(second()-t1-tm2)/ntimes;
dgesl(aa,ldaa,n,ipvt,b,0);
time[1][7]=(second()-t1)/ntimes;
total=time[0][7]+time[1][7];
time[2][7]=total;
time[3][7]=ops/(1.0e3*total);
time[4][7]=2.0e3/time[3][7];
time[5][7]=total/cray;
/*thefollowingcodesequenceimplementsthesemanticsof
theFortranintrinsics"
nint(min(time[3][3],time[3][7]))"
*/
kf=(time[3][3]<
time[3][7])?
time[3][3]:
time[3][7];
kf=(kf>
ZERO)?
(kf+.5):
(kf-.5);
if(fabs((double)kf)<
ONE)
kflops=0;
else{
kflops=floor(fabs((double)kf));
if(kf<
ZERO)kflops=-kflops;
timesforarraywithleadingdimensionof%4d\n"
ldaa);
print_time(4);
print_time(5);
print_time(6);
print_time(7);
Precision%5dKflops;
%dReps\n"
kflops,NTIMES);
}
/*----------------------*/
print_time(row)
introw;
%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n"
(double)time[0][row],
(double)time[1][row],(double)time[2][row],(double)time[3][row],
(double)time[4][row],(double)time[5][row]);
matgen(a,lda,n,b,norma)
REALa[],b[],*norma;
intlda,n;
/*Wewouldliketodeclarea[][lda],butcdoesnotallowit.Inthis
function,referencestoa[i][j]arewrittena[lda*j+i].*/
intinit,i,j;
init=1325;
*norma=0.0;
for(j=0;
j<
j++){
for(i=0;
init=3125*init%65536;
a[lda*j+i]=(init-32768.0)/16384.0;
*norma=(a[lda*j+i]>
*norma)?
a[lda*j+i]:
*norma;
}
b[i]=0.0;
b[i]=b[i]+a[lda*j+i];
dgefa(a,lda,n,ipvt,info)
REALa[];
intlda,n,ipvt[],*info;
function,referencestoa[i][j]arewrittena[lda*i+j].*/
dgefafactorsadoubleprecisionmatrixbygaussianelimination.
dgefaisusuallycalledbydgeco,butitcanbecalled
directlywithasavingintimeifrcondisnotneeded.
(timefordgeco)=(1+9/n)*(timefordgefa).
onentry
aREALprecision[n][lda]
thematrixtobefactored.
ldainteger
theleadingdimensionofthearraya.
ninteger
theorderofthematrixa.
onreturn
aanuppertriangularmatrixandthemultipliers
whichwereusedtoobtainit.
thefactorizationcanbewrittena=l*uwhere
lisaproductofpermutationandunitlower
triangularmatricesanduisuppertriangular.
ipvtinteger[n]
anintegervectorofpivotindices.
infointeger
=0normalvalue.
=kifu[k][k].eq.0.0.thisisnotanerror
conditionforthissubroutine,butitdoes
indicatethatdgeslordgediwilldividebyzero
ifcalled.usercondindgecoforareliable
indicationofsingularity.
linpack.thisversiondated08/14/78.
clevemoler,universityofnewmexico,argonnenationallab.
functions
blasdaxpy,dscal,idamax
*/
/*internalvariables*/
REALt;
intidamax(),j,k,kp1,l,nm1;
/*gaussianeliminationwithpartialpivoting*/
*info=0;
nm1=n-1;
if(nm1>
=0){
for(k=0;
k<
nm1;
k++){
kp1=k+1;
/*findl=pivotindex*/
l=idamax(n-k,&
a[lda*k+k],1)+k;
ipvt[k]=l;
/*zeropivotimpliesthiscolumnalready
triangularized*/
if(a[lda*k+l]!
=ZERO){
/*interchangeifnecessary*/
if(l!
=k){
t=a[lda*k+l];
a[lda*k+l]=a[lda*k+k];
a[lda*k+k]=t;
}
/*computemultipliers*/
t=-ONE/a[lda*k+k];
dscal(n-(k+1),t,&
a[lda*k+k+1],1);
/*roweliminationwithcolumnindexing*/
for(j=kp1;
t=a[lda*j+l];
if(l!
a[lda*j+l]=a[lda*j+k];
a[lda*j+k]=t;
}
daxpy(n-(k+1),t,&
a[lda*k+k+1],1,
&
a[lda*j+k+1],1);
}
else{
*info=k;
}
}
ipvt[n-1]=n-1;
if(a[lda*(n-1)+(n-1)]==ZERO)*info=n-1;
dgesl(a,lda,n,ipvt,b,job)
intlda,n,ipvt[],job;
REALa[],b[];
dgeslsolvesthedoubleprecisionsystem
a*x=bortrans(a)*x=b
usingthefactorscomputedbydgecoordgefa.
adoubleprecision[n][lda]
theoutputfromdgecoordgefa.
thepivotvectorfromdgecoordgefa.
bdoubleprecision[n]
therighthandsidevector.
jobinteger
=0tosolvea*x=b,
=nonzerotosolvetrans(a)*x=bwhere
trans(a)isthetranspose.
bthesolutionvectorx.
errorcondition
adivisionbyzerowilloccuriftheinputfactorcontainsa
zeroonthediagonal.technicallythisindicatessingularity