linpackdoc.docx
《linpackdoc.docx》由会员分享，可在线阅读，更多相关《linpackdoc.docx（25页珍藏版）》请在冰豆网上搜索。
linpackdoc
linpack
/*
TranslatedtoCbyBonnieToy5/88
-modifiedon2/25/94tofixaproblemwithdaxpyfor
unequalincrementsorequalincrementsnotequalto1.
JackDongarra
-modifiedon08/27/09fixtypoline270,plusset'ix'to0
inthecaseincxisnot1
JulieLangou
TocompilesingleprecisionversionforSun-4:
cc-DSP-O4-fsingle-fsingle2clinpack.c-lm
TocompiledoubleprecisionversionforSun-4:
cc-DDP-O4clinpack.c-lm
ToobtainrolledsourceBLAS,add-DROLLto
#defineONE1.0e0
#definePREC"Double"
#endif
#defineNTIMES10
#ifdefROLL
#defineROLLING"Rolled"
#endif
#ifdefUNROLL
#defineROLLING"Unrolled"
#endif
#include
#include
staticREALtime[9][9];
main（）
{
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,"PrecisionLinpack\n\n"）;
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;
t1=second（）;
dgesl（a,lda,n,ipvt,b,0）;
time[1][0]=second（）-t1;
total=time[0][0]+time[1][0];
/*computearesidualtoverifyresults.*/
for（i=0;ix[i]=b[i];
}
matgen（a,lda,n,b,&norma）;
for（i=0;ib[i]=-b[i];
}
dmxpy（n,b,n,lda,x,a）;
resid=0.0;
normx=0.0;
for（i=0;iresid=（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"）;
printf（"%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"）;
fprintf（stderr,"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;
fprintf（stderr,"timesforarraywithleadingdimensionof%5d\n",lda）;
print_time（0）;
matgen（a,lda,n,b,&norma）;
t1=second（）;
dgefa（a,lda,n,ipvt,&info）;
time[0][1]=second（）-t1;
t1=second（）;
dgesl（a,lda,n,ipvt,b,0）;
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;
matgen（a,lda,n,b,&norma）;
t1=second（）;
dgefa（a,lda,n,ipvt,&info）;
time[0][2]=second（）-t1;
t1=second（）;
dgesl（a,lda,n,ipvt,b,0）;
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;
t1=second（）;
for（i=0;itm=second（）;
matgen（a,lda,n,b,&norma）;
tm2=tm2+second（）-tm;
dgefa（a,lda,n,ipvt,&info）;
}
time[0][3]=（second（）-t1-tm2）/ntimes;
t1=second（）;
for（i=0;idgesl（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,&norma）;
t1=second（）;
dgefa（aa,ldaa,n,ipvt,&info）;
time[0][4]=second（）-t1;
t1=second（）;
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;
matgen（aa,ldaa,n,b,&norma）;
t1=second（）;
dgefa（aa,ldaa,n,ipvt,&info）;
time[0][5]=second（）-t1;
t1=second（）;
dgesl（aa,ldaa,n,ipvt,b,0）;
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;
matgen（aa,ldaa,n,b,&norma）;
t1=second（）;
dgefa（aa,ldaa,n,ipvt,&info）;
time[0][6]=second（）-t1;
t1=second（）;
dgesl（aa,ldaa,n,ipvt,b,0）;
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;
ntimes=NTIMES;
tm2=0;
t1=second（）;
for（i=0;itm=second（）;
matgen（aa,ldaa,n,b,&norma）;
tm2=tm2+second（）-tm;
dgefa（aa,ldaa,n,ipvt,&info）;
}
time[0][7]=（second（）-t1-tm2）/ntimes;
t1=second（）;
for（i=0;idgesl（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][3]:
time[3][7];
kf=（kf>ZERO）?
（kf+.5）:
（kf-.5）;
if（fabs（（double）kf）kflops=0;
else{
kflops=floor（fabs（（double）kf））;
if（kf}
fprintf（stderr,"timesforarraywithleadingdimensionof%4d\n",ldaa）;
print_time（4）;
print_time（5）;
print_time（6）;
print_time（7）;
fprintf（stderr,ROLLING）;fprintf（stderr,PREC）;
fprintf（stderr,"Precision%5dKflops;%dReps\n",kflops,NTIMES）;
}
/*----------------------*/
print_time（row）
introw;
{
fprintf（stderr,"%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;jfor（i=0;iinit=3125*init%65536;
a[lda*j+i]=（init-32768.0）/16384.0;
*norma=（a[lda*j+i]>*norma）?
a[lda*j+i]:
*norma;
}
}
for（i=0;ib[i]=0.0;
}
for（j=0;jfor（i=0;ib[i]=b[i]+a[lda*j+i];
}
}
}
/*----------------------*/
dgefa（a,lda,n,ipvt,info）
REALa[];
intlda,n,ipvt[],*info;
/*Wewouldliketodeclarea[][lda],butcdoesnotallowit.Inthis
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;kkp1=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;jt=a[lda*j+l];
if（l!
=k）{
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[];
/*Wewouldliketodeclarea[][lda],butcdoesnotallowit.Inthis
function,referencestoa[i][j]arewrittena[lda*i+j].*/
/*
dgeslsolvesthedoubleprecisionsystem
a*x=bortrans（a）*x=b
usingthefactorscomputedbydgecoordgefa.
onentry
adoubleprecision[n][lda]
theoutputfromdgecoordgefa.
ldainteger
theleadingdimensionofthearraya.
ninteger
theorderofthematrixa.
ipvtinteger[n]
thepivotvectorfromdgecoordgefa.
bdoubleprecision[n]
therighthandsidevector.
jobinteger
=0tosolvea*x=b,
=nonzerotosolvetrans（a）*x=bwhere
trans（a）isthetranspose.
onreturn
bthesolutionvectorx.
errorcondition
adivisionbyzerowilloccuriftheinputfactorcontainsa
zeroonthediagonal.technicallythisindicatessingularity
升级会员 