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
升级会员 