C-------------------------------------------------------------------------- C THE CODE IS WRITTEN BASED ON THE MLPG FORMULATION FOR THE ELASTO- C DYNAMIC ANALYSIS C C MODIFIED FROM NCETEN1.F FOR INCLUDING THE CENTRAL DIFFERENCE METHOD C MODIFIED FROM NCENAD1.F FOR COMMENTS OF THE CODE C 2534 NODES C NO DIFFRACTION CRITERIA C USING NEWMARK FAMILY OF METHODS C USING THE PENALTY METHOD TO ENFORCE THE ESSENTAIL BOUNDARY CONDITIONS C USING NEAR-TIP STRESS FIELDS TO COMPUTE STRESS INTENSITY FACTORS C USING NEWMARK FAMILY OF METHODS FOR TIME INTEGRATION C C USING LAPACK & NAPACK SUBROUTINES C-------------------------------------------------------------------------- C PROGRAM MESHLESS C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /ITYPE/ IOPT,IWEI COMMON /INTEG/ ITIME COMMON /COF/ EMODULE,POISSON COMMON /MBO/ M COMMON /SCA/ SCALE COMMON /FWE/ KW COMMON /GAU1/ NGL,NGLX,NGLY COMMON /PEN/ PENA COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /SPEC/ NDSPX(2534),NDSPY(2534),NLDX(2534),NLDY(2534) COMMON /NIPT/ SDSPX(2534),SDSPY(2534),SLDX(2534),SLDY(2534) COMMON /BOUND/ NBOUND,IBOUND(2534),INDBND(2534) COMMON /DEN/ CT COMMON /PA/ ALFA,GAMA COMMON /NST/ NTIME,STIME,FTIME COMMON /TS/ DT,A1,A2,A3,A4,A5 COMMON /CRACK/ CX1,CY1,CX2,CY2 COMMON /NOUP/ NRSE,NSG(6) COMMON /OUTT/ NOUPT COMMON /OUPD/ SL,NINP COMMON /OUTK/ NEX COMMON /RANGE/ XK1A,XK1B,XK2A,XK2B COMMON /LOAD/ XLOAD COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) COMMON /CON/ D(3,3) COMMON /GAU2/ POINT2(9,2),WEIGHT2(9,2) COMMON /ROCCR/ RO(2534),CC(2534),R(2534) COMMON /INTANG/ THETA1(2534),THETA2(2534) C DIMENSION SKM(5068,5068),FM(5068),FSOL(5068) DIMENSION SK(5068,5068),SM(5068),F(5068) DIMENSION POINT1(9),WEIGHT1(9) C CALL INPUT1 CALL FEMATIS CALL FIN2 CALL GRO CALL ANGLE C DO 50 IST=1,NTIME+1 C ICT=IST-1 STIST=ICT*DT C WRITE(*,*) ICT,STIST C CALL INSKF(SKM,FM,FSOL) C IF (IST .EQ. 1) THEN C CALL INSKF(SK,SM,F) C DO 10 IEL=1,NNODE C WRITE(*,*) IEL C CALL XYQPO(SK,SM,IEL) C 10 CONTINUE C NGLT=NGL CALL FIN1(POINT1,WEIGHT1,NGLT) CALL BCAPLY(SK,F,POINT1,WEIGHT1) C CALL IACCE(SM,F) C IF (ITIME .EQ. 1) THEN C CALL TIMESKM(SKM,SK,SM) C CALL INVEKM(SKM,SK) C END IF C ELSE C IF (ITIME .EQ. 1) THEN C CALL TIMEFM(FM,F,SM) C CALL SOLVE(SK,FM,FSOL) C CALL NVEAC(FSOL) C ELSE IF (ITIME .EQ. 2) THEN C CALL TIMER(FM,FSOL,SK,SM,F) C END IF C C FOR THE OUTPUT OF THE RESULTS C CALL OUTP(ICT,STIST) C IF (ICT .EQ. NOUPT) THEN C CALL OUTALL(ICT,STIST) CALL OUTPA(ICT,STIST) CALL DEFO(ICT,STIST) C END IF C C COMPUTE STRESS INTENSITY FACTORS BY EXTRAPOLATION C CALL KCOMP(ICT,STIST) END IF C 50 CONTINUE C STOP C END C C C---------------------------------------------------------------------------------- C INPUT1 C C DESCRIPTION OF THE VARIABLES C C IOPT = 1:PLANE STRESS ; IOPT = 2:PLANE STRAIN C IEWI = 1:GAUSSIAN WEIGHT FUNCTION ; IEWI = 2:SPLINE WEIGHT FUNCTION C ITIME = 1: NEWMARK METHOD, ITIME = 2 CENTRAL DIFFERENCE METHOD C EMODULE,POISSON : YOUNG'S MODULUS, POISSON'S RATIO (RESPECTIVELY) C M : ORDER OF THE DISPLACEMENT BASIS VECTOR C M = 3:LINEAR BASIS ; M = 6:QUADRATIC BASIS C SCALE : DILATION (SCALING FACTOR) FOR DETERMINING THE DOMAIN OF INFLUENCE C KW : THE POWER OF THE EXPONENT IN THE GAUSSIAN WEIGHT FUNCTION C NGL : NUMBER OF GAUSSIAN INTEGRATION POINTS FOR THE BOUNDARY INTEGRAL C NGLX,NGLY : NUMBER OF X,Y (RESCEPTIVELY) GAUSSIAN INTEGRATION POINTS C FOR THE DOMAIN INTEGRAL C PENA : PENALTY PARAMETER USED TO ENFORCE ESSENTIAL BCs C NNODE : TOTAL NUMBER OF NODES C X(I),Y(I) : X,Y (RESPECTIVELY) POSITIONS C NDSPX(I),NDSPY(I) : NODAL SPECIFIED DISPLACEMENT FLAG C >0 NON SPECIFIED (FREE), <0 SPECIFID (=-2 FIRST NODE C SPECIFIED, =-3 LAST NODE SPECIFIED) (CCW) C NLDX(I),NLDY(I) : NODAL SPECIFIED LOAD FLAG C >0 NON SPECIFIED (FREE), <0 SPECIFID (=-2 FIRST NODE C SPECIFIED, =-3 LAST NODE SPECIFIED) (CCW) C SDSPX(I),SDSPY(I) : SPECIFIED X,Y DISPLACEMENTS (RESPECTIVELY) C SLDX(I),SLDY(I) : SPECIFIED X,Y LOADS (RESPECTIVELY) C NBOUND : TOTAL NUMBER OF BOUNDARY NODES C INDBND(I) : THE NODE NUMBER FOR THE Ith NODE ALONG THE BOUNDARY C IBOUND(I*) : THE Ith NODE CORRESPONDING TO THE NODE NUMBER I* C C *NOTE : BOUNDARY NODES MUST BE ENTERED IN A COUNTER-CLOCKWISE C ORDER SO THAT THE NORMALS AND ANGLES ARE PROPERLY CALCULATED C C CT : THE MASS DENSITY C ALFA,GAMA : PARAMETERS IN THE NEWMARK SCHEME C ALFA=1/2, GAMA=1/2 (CONSTANT AVERAGE ACCELERATION METHOD) C ALFA=1/2, GAMA=1/3 (LINEAR ACCELERATION METHOD) C ALFA=3/2, GAMA=8/5 (GALERKIN METHOD) C ALFA=3/2, GAMA=2 (BACKWARD DIFFERENCE METHOD) C NTIME, STIME, FTIME : NTIME = NUMBER OF TIME STEPS C STIME = THE STARTING TIME C FTIME = THE ENDING TIME C DT : THE TIME STEP C CX1,CY1,CX2,CY2 : THE X,Y COORDIANTES OF BOTH END POINTS OF THE CRACK LINE C CX1,CY1 = LEFT END POINT C CX2,CY2 = RIGHT END POINT C NRSE,NSG(I) : NUMBER OF THE NODES DIRECTLY AHEAD OF THE CRACK TIP AND THEIR C NODE NO. (RESPECTIVELY) C NOUPT : THE SPECIFIED TIME STEP FOR OUTPUT C SL, NINP : RADIAL DISTANCE ALONG THE CRACK LINE AND NUMBER OF POINTS C FOR OUTPUT (RESPECTIVELY) C NEX : NUMBER OF POINTS FOR COMPUTING STRESS INTENSITY FACTOR ALONG A RANGE C (R) NEAR THE CRACK TIP C XK1A,XK1B,XK2A,XK2B : XK1A,XK1B = THE STARTING AND ENDING POINTS FOR THE C FIRST R RANGE C XK2A,XK2B = THE STARTING AND ENDING POINTS FOR THE C SECOND R RANGE C XLOAD : THE TRACTION FORCES (ASSUMED TO BE UNIFORM) C----------------------------------------------------------------------------------- C SUBROUTINE INPUT1 C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*81 DUMMY C COMMON /ITYPE/ IOPT,IWEI COMMON /INTEG/ ITIME COMMON /COF/ EMODULE,POISSON COMMON /MBO/ M COMMON /SCA/ SCALE COMMON /FWE/ KW COMMON /GAU1/ NGL,NGLX,NGLY COMMON /PEN/ PENA COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /SPEC/ NDSPX(2534),NDSPY(2534),NLDX(2534),NLDY(2534) COMMON /NIPT/ SDSPX(2534),SDSPY(2534),SLDX(2534),SLDY(2534) COMMON /BOUND/ NBOUND,IBOUND(2534),INDBND(2534) COMMON /DEN/ CT COMMON /PA/ ALFA,GAMA COMMON /NST/ NTIME,STIME,FTIME COMMON /TS/ DT,A1,A2,A3,A4,A5 COMMON /CRACK/ CX1,CY1,CX2,CY2 COMMON /NOUP/ NRSE,NSG(6) COMMON /OUTT/ NOUPT COMMON /OUPD/ SL,NINP COMMON /OUTK/ NEX COMMON /RANGE/ XK1A,XK1B,XK2A,XK2B COMMON /LOAD/ XLOAD COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C LOGICAL ISDSP,ISLD C OPEN (UNIT=7,FILE='ncein4.dat',STATUS='OLD') OPEN (UNIT=8,FILE='in-chec.dat') C READ (7,'(A)') DUMMY READ (7,'(A)') DUMMY READ (7,*) IOPT,IWEI WRITE (8,*) IOPT,IWEI C READ (7,'(A)') DUMMY READ (7,*) ITIME WRITE (8,*) ITIME C READ (7,'(A)') DUMMY READ (7,*) EMODULE,POISSON WRITE (8,*) EMODULE,POISSON C IF (IOPT .EQ. 2) THEN EMODULE=EMODULE/(1.0-POISSON**2) POISSON=POISSON/(1.0-POISSON) END IF C WRITE (8,*) EMODULE,POISSON C READ (7,'(A)') DUMMY READ (7,*) M WRITE (8,*) M C READ (7,'(A)') DUMMY READ (7,*) SCALE WRITE (8,*) SCALE C READ (7,'(A)') DUMMY READ (7,*) KW WRITE (8,*) KW C READ (7,'(A)') DUMMY READ (7,*) NGL,NGLX,NGLY WRITE (8,*) NGL,NGLX,NGLY C READ (7,'(A)') DUMMY READ (7,*) PENA WRITE (8,*) PENA C READ (7,'(A)') DUMMY READ (7,*) NNODE DO 10 I=1,NNODE READ (7,*) IC,X(I),Y(I),NDSPX(I),NDSPY(I),NLDX(I),NLDY(I) 10 CONTINUE C WRITE (8,*) NNODE DO 15 I=1,NNODE WRITE (8,*) I,X(I),Y(I),NDSPX(I),NDSPY(I),NLDX(I),NLDY(I) 15 CONTINUE C READ (7,'(A)') DUMMY DO 20 I=1,NNODE C ISDSP=NDSPX(I).LT.0.OR.NDSPY(I).LT.0 ISLD=NLDX(I).LT.0.OR.NLDY(I).LT.0 C IF (ISDSP .OR. ISLD) THEN READ (7,*) IC,SDSPX(IC),SDSPY(IC),SLDX(IC),SLDY(IC) END IF C 20 CONTINUE C DO 30 I=1,NNODE INDBND(I)=0 IBOUND(I)=0 30 CONTINUE C READ (7,'(A)') DUMMY READ (7,*) NBOUND DO 40 I=1,NBOUND READ (7,*) INDBND(I) IBOUND(INDBND(I))=I 40 CONTINUE C C READ THE NECESSARY DATA FOR TIME-DEPENDENT PROBLEMS C READ (7,'(A)') DUMMY READ (7,*) CT WRITE (8,*) CT C READ (7,'(A)') DUMMY READ (7,*) ALFA,GAMA WRITE (8,*) ALFA,GAMA C READ (7,'(A)') DUMMY READ (7,*) NTIME,STIME,FTIME WRITE (8,*) NTIME,STIME,FTIME C DT=(FTIME-STIME)/NTIME DT2=DT*DT C A1=ALFA*DT A2=(1.0-ALFA)*DT A3=2.0/GAMA/DT2 A4=A3*DT A5=1.0/GAMA-1.0 C WRITE (8,*) DT,A1,A2,A3,A4,A5 C C END OF DATA FOR TIME-DEPENDENT PROBLEMS C READ (7,'(A)') DUMMY READ (7,*) CX1,CY1,CX2,CY2 WRITE (8,*) CX1,CY1,CX2,CY2 C READ (7,'(A)') DUMMY READ (7,*) NRSE WRITE(8,*) NRSE C DO 50 I=1,NRSE READ (7,*) NSG(I) 50 CONTINUE DO 55 I=1,NRSE WRITE (8,*) NSG(I) 55 CONTINUE C READ (7,'(A)') DUMMY READ (7,*) NOUPT WRITE (8,*) NOUPT C READ (7,'(A)') DUMMY READ (7,*) SL,NINP WRITE (8,*) SL,NINP C READ (7,'(A)') DUMMY READ (7,*) NEX READ (7,*) XK1A,XK1B,XK2A,XK2B READ (7,*) XLOAD C WRITE (8,*) NEX WRITE (8,*) XK1A,XK1B,XK2A,XK2B WRITE (8,*) XLOAD C C INITIALIZE THE GLOBAL DISPLACEMENT, VELOCITY AND ACCELERATION C DO 60 I=1,2*NNODE GLU(I)=0.0 GLV(I)=0.0 GLA(I)=0.0 60 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C FEMATIS C C CREATE THE 3X3 CONSTITUTIVE MATRIX (D) FOR AN ISOTROPIC C LINEARLY ELASTIC MATERIAL C------------------------------------------------------------------------ C SUBROUTINE FEMATIS C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /COF/ EMODULE,POISSON COMMON /CON/ D(3,3) C CONST=EMODULE/(1.0-POISSON*POISSON) D(1,1)=CONST D(1,2)=CONST*POISSON D(1,3)=0.0 D(2,1)=D(1,2) D(2,2)=D(1,1) D(2,3)=0.0 D(3,1)=0.0 D(3,2)=0.0 D(3,3)=CONST*(1.0-POISSON)/2.0 C RETURN END C C C------------------------------------------------------------------------ C FIN2 C C LOCATE THE GAUSSIAN INTEGRATION POINTS IN THE NATURAL PLANE AND C ASSIGN WEIGHTS OVER A SQUARE DOMAIN (-1<=SKI<=1, -1<=ETA<=1) C------------------------------------------------------------------------ C SUBROUTINE FIN2 C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /GAU1/NGL,NGLX,NGLY COMMON /GAU2/POINT2(9,2),WEIGHT2(9,2) C DIMENSION POINTX(9),POINTY(9),WEIGHTX(9),WEIGHTY(9) C IF (NGLX .GT. NGLY) THEN NGLI=NGLX ELSE NGLI=NGLY END IF C DO 10 I=1,NGLI DO 10 J=1,2 POINT2(I,J)=0.0 WEIGHT2(I,J)=0.0 10 CONTINUE C NGLT=NGLX CALL FIN1(POINTX,WEIGHTX,NGLT) C NGLT=NGLY CALL FIN1(POINTY,WEIGHTY,NGLT) C DO 20 INTX=1,NGLX POINT2(INTX,1)=POINTX(INTX) WEIGHT2(INTX,1)=WEIGHTX(INTX) 20 CONTINUE C DO 30 INTY=1,NGLY POINT2(INTY,2)=POINTY(INTY) WEIGHT2(INTY,2)=WEIGHTY(INTY) 30 CONTINUE C RETURN END C C C------------------------------------------------------------------------- C FIN1 C C LOCATE THE 1D GAUSSIAN INTEGRATION POINTS IN THE NATURAL PLANE AND C ASSIGN WEIGHTS (SKI) C------------------------------------------------------------------------- C SUBROUTINE FIN1(POINT1,WEIGHT1,NGLT) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION POINT1(9),WEIGHT1(9) C DO 100 I=1,NGLT POINT1(I)=0.0 WEIGHT1(I)=0.0 100 CONTINUE C C FINDING CORRESPONDING INTEGRATION POINTS AND WEIGHTS C IF (NGLT .EQ. 1) THEN POINT1(1)=0.0 WEIGHT1(1)=2.0 C ELSE IF (NGLT .EQ. 2) THEN POINT1(1)=-0.577350269189626 POINT1(2)=-POINT1(1) WEIGHT1(1)=1.0 WEIGHT1(2)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 3) THEN POINT1(1)=-0.774596669241483 POINT1(2)=0.0 POINT1(3)=-POINT1(1) WEIGHT1(1)=0.555555555555556 WEIGHT1(2)=0.888888888888889 WEIGHT1(3)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 4) THEN POINT1(1)=-0.861136311594053 POINT1(2)=-0.339981043584856 POINT1(3)=-POINT1(2) POINT1(4)=-POINT1(1) WEIGHT1(1)=0.347854845137454 WEIGHT1(2)=0.652145154862546 WEIGHT1(3)=WEIGHT1(2) WEIGHT1(4)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 5) THEN POINT1(1)=-0.906179845938664 POINT1(2)=-0.538469310105683 POINT1(3)=0.0 POINT1(4)=-POINT1(2) POINT1(5)=-POINT1(1) WEIGHT1(1)=0.236926885056189 WEIGHT1(2)=0.478628670499366 WEIGHT1(3)=0.568888888888889 WEIGHT1(4)=WEIGHT1(2) WEIGHT1(5)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 6) THEN POINT1(1)=0.932469514203252 POINT1(2)=0.661209386466265 POINT1(3)=0.238619186083197 POINT1(4)=-POINT1(3) POINT1(5)=-POINT1(2) POINT1(6)=-POINT1(1) WEIGHT1(1)=0.171324492379107 WEIGHT1(2)=0.360761573048139 WEIGHT1(3)=0.467913934572691 WEIGHT1(4)=WEIGHT1(3) WEIGHT1(5)=WEIGHT1(2) WEIGHT1(6)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 7) THEN POINT1(1)=0.949107912342759 POINT1(2)=0.741531185599394 POINT1(3)=0.405845151377397 POINT1(4)=0.0 POINT1(5)=-POINT1(3) POINT1(6)=-POINT1(2) POINT1(7)=-POINT1(1) WEIGHT1(1)=0.129484966168870 WEIGHT1(2)=0.279705391489277 WEIGHT1(3)=0.381830050505119 WEIGHT1(4)=0.417949183673469 WEIGHT1(5)=WEIGHT1(3) WEIGHT1(6)=WEIGHT1(2) WEIGHT1(7)=WEIGHT1(1) C ELSE IF (NGLT .EQ. 8) THEN POINT1(1)=0.960289856497536 POINT1(2)=0.796666477413627 POINT1(3)=0.525532409916329 POINT1(4)=0.183434642495650 POINT1(5)=-POINT1(4) POINT1(6)=-POINT1(3) POINT1(7)=-POINT1(2) POINT1(8)=-POINT1(1) WEIGHT1(1)=0.101228536290376 WEIGHT1(2)=0.222381034423374 WEIGHT1(3)=0.313706645877887 WEIGHT1(4)=0.362683783378362 WEIGHT1(5)=WEIGHT1(4) WEIGHT1(6)=WEIGHT1(3) WEIGHT1(7)=WEIGHT1(2) WEIGHT1(8)=WEIGHT1(1) C ELSE POINT1(1)=0.968160239507626 POINT1(2)=0.836061107326636 POINT1(3)=0.613371432700590 POINT1(4)=0.324253423403809 POINT1(5)=0.0 POINT1(6)=-POINT1(4) POINT1(7)=-POINT1(3) POINT1(8)=-POINT1(2) POINT1(9)=-POINT1(1) WEIGHT1(1)=0.081274388361574 WEIGHT1(2)=0.180648160694857 WEIGHT1(3)=0.260610696402935 WEIGHT1(4)=0.312347077040003 WEIGHT1(5)=0.330239355001260 WEIGHT1(6)=WEIGHT1(4) WEIGHT1(7)=WEIGHT1(3) WEIGHT1(8)=WEIGHT1(2) WEIGHT1(9)=WEIGHT1(1) C END IF C RETURN END C C C---------------------------------------------------------------------------- C GRO C C CALCULATE THE DOMAIN OF INFLUENCE AND DOMAIN OF INTEGRATION FOR C EVERY NODE AND FINDS THE DISTANCE TO THE NEAREST NEIGHBOR NODE C C RO(I) : THE RADIUS OF THE LOCAL SUBDOMAIN C R(I) : THE RADIUS OF THE DOMAIN OF INFLUENCE C CC(I) : THE DISTANCE TO THE THIRD NEAREST NEIGHBORING POINT FROM C NODE I C SCALE : DILATION (SCALING FACTOR) FOR DETERMING THE DOMAIN OF INFLUENCE C R=SCALE*CC C---------------------------------------------------------------------------- C SUBROUTINE GRO C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /SCA/ SCALE COMMON /ROCCR/ RO(2534),CC(2534),R(2534) C REAL*4 SORTD(NNODE-1) DIMENSION DIS(NNODE-1),DW(NNODE-1) C DO 10 IEL=1,NNODE RO(IEL)=0.0 CC(IEL)=0.0 R(IEL)=0.0 10 CONTINUE C NNODE1=NNODE-1 C DO 20 IEL=1,NNODE C KK=0 DO 30 IDL=1,NNODE1 DIS(IDL)=0.0 30 CONTINUE C DO 40 IN=1,NNODE IF (IN .NE. IEL) THEN KK=KK+1 XD=X(IN)-X(IEL) YD=Y(IN)-Y(IEL) DIS(KK)=DSQRT(XD**2+YD**2) END IF 40 CONTINUE C C CALL NAPACK SUBROUTINES C CALL SORT2(DIS,SORTD,DW,NNODE1) C C DIST1=DIS(SORTD(1)) C DIST2=DIS(SORTD(2)) C J=2 C C DO WHILE (DIST1 .EQ. DIST2) C J=J+1 C DIST2=DIS(SORTD(J)) C END DO C C DIST3=DIS(SORTD(J+1)) C C DO WHILE (DIST2 .EQ. DIST3) C J=J+1 C DIST3=DIS(SORTD(J)) C END DO C C WRITE(12,*) DIST1,DIST2,DIST3 C C RO(IEL)=DIST1 C CC(IEL)=DIST3 C R(IEL)=SCALE*CC(IEL) C RO(IEL)=DIS(SORTD(1)) CC(IEL)=DIS(SORTD(3)) R(IEL)=SCALE*CC(IEL) C 20 CONTINUE C RETURN END C C C----------------------------------------------------------------------- C ANGLE C C CALCULATE ANGLES (THETA1,THETA2) FOR THE NODAL DOMAIN OF C INTEGRATION. FOR INTERNAL NODES (I.E., NOT ON BOUNDARY) THESE ARE C ZERO AND 2PI. FOR ALL OTHER NODES (BOUNDARIES) THE START AND C FINISH ANGLES ARE DETERMINED BASED UPON THE GEOMETRY ASSUMING A C LINEAR BOUNDARY CURVE FIT BETWEEN NODES. C----------------------------------------------------------------------- C SUBROUTINE ANGLE C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /NIPT/ SDSPX(2534),SDSPY(2534),SLDX(2534),SLDY(2534) COMMON /BOUND/ NBOUND,IBOUND(2534),INDBND(2534) COMMON /INTANG/ THETA1(2534),THETA2(2534) C DATA PI /3.141592653/ C DO 10 I=1,NNODE THETA1(I)=0.0 THETA2(I)=2*PI 10 CONTINUE C DO 100 I=1,NBOUND C C FINDING VECTOR TO NEXT (I+1) AND PREVIOUS (I-1) BOUNDARY NODES C IF (I .EQ. 1) THEN R1X = X(INDBND(I+1))-X(INDBND(I)) R1Y = Y(INDBND(I+1))-Y(INDBND(I)) R2X = X(INDBND(NBOUND))-X(INDBND(I)) R2Y = Y(INDBND(NBOUND))-Y(INDBND(I)) ELSE IF (I .EQ. NBOUND) THEN R1X = X(INDBND(1))-X(INDBND(I)) R1Y = Y(INDBND(1))-Y(INDBND(I)) R2X = X(INDBND(I-1))-X(INDBND(I)) R2Y = Y(INDBND(I-1))-Y(INDBND(I)) ELSE R1X = X(INDBND(I+1))-X(INDBND(I)) R1Y = Y(INDBND(I+1))-Y(INDBND(I)) R2X = X(INDBND(I-1))-X(INDBND(I)) R2Y = Y(INDBND(I-1))-Y(INDBND(I)) END IF C R1=DSQRT(R1X**2+R1Y**2) R2=DSQRT(R2X**2+R2Y**2) C C CALCULATING BASE ANGLES C THETAA=DACOS(R1X/R1) THETAB=DACOS(R2X/R2) C C FINDING QUADRANT FOR BOTH NODES C IF (R1X.GE.0.0 .AND. R1Y.GE.0.0) THEN NQ1 = 1 ELSE IF (R1X.LE.0.0 .AND. R1Y.GT.0.0) THEN NQ1 = 2 ELSE IF (R1X.LT.0.0 .AND. R1Y.LE.0.0) THEN NQ1 = 3 ELSE NQ1 = 4 END IF C IF (R2X.GE.0.0 .AND. R2Y.GE.0.0) THEN NQ2 = 1 ELSE IF (R2X.LE.0.0 .AND. R2Y.GT.0.0) THEN NQ2 = 2 ELSE IF (R2X.LT.0.0 .AND. R2Y.LE.0.0) THEN NQ2 = 3 ELSE NQ2 = 4 END IF C C CALCULATING TRUE ANGLES C C THETA1 C IF (NQ1.EQ.1 .OR. NQ1.EQ.2) THEN THETA1(INDBND(I))=THETAA ELSE THETA1(INDBND(I))=2*PI-THETAA END IF C C THETA2 C IF ((NQ1.EQ.NQ2).AND.(NQ1.EQ.1.OR.NQ1.EQ.2)) THEN C IF (THETAA.GT.THETAB) THEN THETA2(INDBND(I))=THETAB+2*PI ELSE THETA2(INDBND(I))=THETAB END IF C ELSE IF ((NQ1.EQ.NQ2).AND.(NQ1.EQ.3.OR.NQ1.EQ.4)) THEN C IF (THETAA.GT.THETAB) THEN THETA2(INDBND(I))=2*PI-THETAB ELSE THETA2(INDBND(I))=4*PI-THETAB END IF C ELSE IF (NQ1.EQ.1.AND.NQ2.EQ.2) THEN THETA2(INDBND(I))=THETAB ELSE IF (NQ1.EQ.1.AND.(NQ2.EQ.3.OR.NQ2.EQ.4)) THEN THETA2(INDBND(I))=2*PI-THETAB ELSE IF (NQ1.EQ.2.AND.NQ2.EQ.1) THEN THETA2(INDBND(I))=THETAB+2*PI ELSE IF (NQ1.EQ.2.AND.(NQ2.EQ.3.OR.NQ2.EQ.4)) THEN THETA2(INDBND(I))=2*PI-THETAB ELSE IF (NQ1.EQ.3.AND.(NQ2.EQ.1.OR.NQ2.EQ.2)) THEN THETA2(INDBND(I))=THETAB+2*PI ELSE IF (NQ1.EQ.3.AND.NQ2.EQ.4) THEN THETA2(INDBND(I))=2*PI-THETAB ELSE IF (NQ1.EQ.4.AND.(NQ2.EQ.1.OR.NQ2.EQ.2)) THEN THETA2(INDBND(I))=THETAB+2*PI ELSE IF (NQ1.EQ.4.AND.NQ2.EQ.3) THEN THETA2(INDBND(I))=4*PI-THETAB END IF C 100 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C INSKF C C INITIALIZE THE STIFFNESS MATRIX (SK), THE MASS MATRIX (SM), AND THE C LOADING VECTOR (F) C C *NOTE : SM IS DIAGONALIZED BY THE ROW-SUM SCHEME C------------------------------------------------------------------------ C SUBROUTINE INSKF(SK,SM,F) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE C DIMENSION SK(2*NNODE,2*NNODE),SM(2*NNODE),F(2*NNODE) C NEQ=2*NNODE C DO 10 I=1,NEQ F(I)=0.0 SM(I)=0.0 DO 20 J=1,NEQ SK(I,J)=0.0 20 CONTINUE 10 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C XYQPO C C CALCULATE THE DOMAIN PORTION OF STIFFNESS (SK) AND MASS (SM) C------------------------------------------------------------------------ C SUBROUTINE XYQPO(SK,SM,IEL) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /GAU1/ NGL,NGLX,NGLY COMMON /NOTN/ NNODE COMMON /GAU2/POINT2(9,2),WEIGHT2(9,2) COMMON /ROCCR/ RO(2534),CC(2534),R(2534) COMMON /GIPT/ X(2534),Y(2534) C DIMENSION SK(2*NNODE,2*NNODE),SM(2*NNODE) DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION STREST(2,3) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) C ROI=RO(IEL) C DO 10 INTX=1,NGLX C SKI=POINT2(INTX,1) WTX=WEIGHT2(INTX,1) C DO 20 INTY=1,NGLY C ETA=POINT2(INTY,2) WTY=WEIGHT2(INTY,2) C CALL XYQDJ(XQ,YQ,DGJ,IEL,ROI,SKI,ETA) C CALL WCOUNT(XQ,YQ,NN,ND,W,DWX,DWY) C CALL TESFUN(XQ,YQ,IEL,VTI,STREST) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XQ,YQ) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL NSTIFF(SK,SM,IEL,NN,ND,VTI,STREST,PHI,DPHIX,DPHIY,WTX & ,WTY,DGJ) C 20 CONTINUE 10 CONTINUE C RETURN END C C C-------------------------------------------------------------------------------- C XYQDJ C C CALCULATE THE GLOBAL COORDINATE (XQ,YQ) OF THE GIVEN GAUSS POINT (SKI,ETA). C ALSO DETERMINE THE ABSOLUTE VALUE OF THE DETERMINANT OF THE JACOBIAN C (DGJ) C-------------------------------------------------------------------------------- C SUBROUTINE XYQDJ(XQ,YQ,DGJ,IEL,ROI,SKI,ETA) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /GIPT/ X(2534),Y(2534) COMMON /INTANG/ THETA1(2534),THETA2(2534) C RADI=(ROI/2.0)*SKI+ROI/2.0 THETA=0.5*THETA1(IEL)*(1-ETA)+0.5*THETA2(IEL)*(1+ETA) DGJ=RADI*(0.5*ROI)*(0.5*(THETA2(IEL)-THETA1(IEL))) C XQ=X(IEL)+RADI*DCOS(THETA) YQ=Y(IEL)+RADI*DSIN(THETA) C RETURN END C C C------------------------------------------------------------------------------ C WCOUNT C C CALCULATE NUMBER OF NODES WHOSE WEIGHT FUNCTIONS ARE LARGE THAN ZERO C AT (XQ,YQ), AND COMPUTE SPATIAL DERIVATIVES OF THESE WEIGHT FUNCTIONS C C NN : NUMBER OF THE NODES WITH THE INFLUENCE AT THE GAUSSIAN POINT (XQ,YQ) C ND(I) : THE NODE NUMBER FOR THE WEIGHT FUNCTION WITH THE INFLUENCE AT C THE GAUSSIAN POINT (XQ,YQ) C W(I),DWX(I),DWY(I) : THE WEIGHT FUNCTION, THE X-DERIVATIVE OF THE WEIGHT C FUNCTION, THE Y-DERIVATIVE OF THE WEIGHT FUNCTION C------------------------------------------------------------------------------ C SUBROUTINE WCOUNT(XQ,YQ,NN,ND,W,DWX,DWY) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /ITYPE/ IOPT,IWEI COMMON /FWE/ KW COMMON /NOTN/ NNODE COMMON /ROCCR/ RO(2534),CC(2534),R(2534) COMMON /GIPT/ X(2534),Y(2534) C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION XD(NNODE),YD(NNODE),D(NNODE),WW(NNODE) C IF (IWEI .EQ. 1) THEN C NN=0 K2=2*KW DO 20 I=1,NNODE RCI=R(I)/CC(I) DEWW=1-DEXP(-(RCI)**K2) XD(I)=XQ-X(I) YD(I)=YQ-Y(I) D(I)=DSQRT(XD(I)**2+YD(I)**2) DCI=D(I)/CC(I) C WW(I)=(DEXP(-(DCI)**K2)-DEXP(-(RCI)**K2))/DEWW C IF (WW(I) .GT. 0) THEN NN=NN+1 W(NN)=WW(I) DWX(NN)=(DEXP(-(DCI)**K2)*(-K2)*(D(I)**(K2-2))*XD(I) & /(CC(I)**K2))/DEWW DWY(NN)=(DEXP(-(DCI)**K2)*(-K2)*(D(I)**(K2-2))*YD(I) & /(CC(I)**K2))/DEWW ND(NN)=I END IF C 20 CONTINUE C ELSE IF (IWEI .EQ. 2) THEN C NN=0 DO 30 I=1,NNODE XD(I)=XQ-X(I) YD(I)=YQ-Y(I) D(I)=DSQRT(XD(I)**2+YD(I)**2) C WW(I)=1-6*(D(I)/R(I))**2+8*(D(I)/R(I))**3-3*(D(I)/R(I))**4 C IF (WW(I) .GT. 0) THEN NN=NN+1 W(NN)=WW(I) DWX(NN)=-12*XD(I)/R(I)**2+24*D(I)*XD(I)/R(I)**3 & -12*D(I)**2*XD(I)/R(I)**4 DWY(NN)=-12*YD(I)/R(I)**2+24*D(I)*YD(I)/R(I)**3 & -12*D(I)**2*YD(I)/R(I)**4 ND(NN)=I END IF C 30 CONTINUE C END IF C RETURN END C C C------------------------------------------------------------------------- C TESFUN C C CALCULATE THE VALUE OF THE TEST FUNCTION (VTI) FOR THE GIVEN C NODE (IEL) AT THE GIVEN GAUSS POINT (XQ,YQ) C C STREST : THE TEST FUNCTION (2X3) DERIVATIVE MATRIX C------------------------------------------------------------------------- C SUBROUTINE TESFUN(XQ,YQ,IEL,VTI,STREST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /ITYPE/ IOPT,IWEI COMMON /FWE/ KW COMMON /ROCCR/ RO(2534),CC(2534),R(2534) COMMON /GIPT/ X(2534),Y(2534) C DIMENSION STREST(2,3) C ROI=RO(IEL) C IF (IWEI .EQ. 2) GOTO 10 C K2=2*KW CCI=CC(IEL) ROC=ROI/CCI DTEWW=1-DEXP(-(ROC)**K2) C XDE=XQ-X(IEL) YDE=YQ-Y(IEL) DE=DSQRT(XDE**2+YDE**2) C DTCI=DE/CCI VTI=(DEXP(-(DTCI)**K2)-DEXP(-(ROC)**K2))/DTEWW DVTXI=(DEXP(-(DTCI)**K2)*(-K2)*(DE**(K2-2))*XDE/(CCI**K2)) & /DTEWW DVTYI=(DEXP(-(DTCI)**K2)*(-K2)*(DE**(K2-2))*YDE/(CCI**K2)) & /DTEWW C GOTO 20 C C 10 XDE=XQ-X(IEL) YDE=YQ-Y(IEL) DE=DSQRT(XDE**2+YDE**2) C VTI=1-6*(DE/ROI)**2+8*(DE/ROI)**3-3*(DE/ROI)**4 DVTXI=-12*XDE/ROI**2+24*DE*XDE/ROI**3-12*DE**2*XDE/ROI**4 DVTYI=-12*YDE/ROI**2+24*DE*YDE/ROI**3-12*DE**2*YDE/ROI**4 C 20 STREST(1,1)=DVTXI STREST(1,2)=0.0 STREST(1,3)=DVTYI STREST(2,1)=0.0 STREST(2,2)=DVTYI STREST(2,3)=DVTXI C RETURN END C C C--------------------------------------------------------------------------------- C INIAB C C INITALIZE THE MATRICES (A,ADX,ADY,B,BDX,BDY) C--------------------------------------------------------------------------------- C SUBROUTINE INIAB(A,B,ADX,ADY,BDX,BDY,NN) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M C DIMENSION A(M,M),B(M,NN),ADX(M,M),ADY(M,M),BDX(M,NN),BDY(M,NN) C DO 10 I=1,M DO 20 J=1,M A(I,J)=0.0 ADX(I,J)=0.0 ADY(I,J)=0.0 20 CONTINUE 10 CONTINUE C DO 30 I=1,M DO 40 J=1,NN B(I,J)=0.0 BDX(I,J)=0.0 BDY(I,J)=0.0 40 CONTINUE 30 CONTINUE C RETURN END C C C----------------------------------------------------------------------------- C ABM C C CALCULATE THE (A) MATRIX, THE (B) MATRIX AND THEIR DERIVATIVES. C C A IS GIVEN BY: C C A=WI(X)P(XI)PT(XI) C C WHERE WI(X) IS THE ITH NODE WEIGHT AT X, P(XI) IS THE BASIS VECTOR C AT THE ITH NODE AND PT(XI) IS ITS TRANSPOSE. THIS IS SUMMED OVER C I WHICH INCLUDES ALL NODES IN THE DOMAIN OF DEFINITION FOR X C C B IS GIVEN BY: C C B=WI(X)P(XI) C C WHERE B IS A SET OF COLUMN VECTORS [WI(X)P(XI)] THAT IS NN WIDE C----------------------------------------------------------------------------- C SUBROUTINE ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M C DIMENSION A(M,M),B(M,NN),ADX(M,M),ADY(M,M),BDX(M,NN),BDY(M,NN) DIMENSION PB2(6),PP2(6,6) DIMENSION ND(NN),W(NN),DWX(NN),DWY(NN) C DO 10 IN=1,NN C CALL BAS(PB2,PP2,IN,ND,NN) C WEI=W(IN) DWEIX=DWX(IN) DWEIY=DWY(IN) C DO 20 I=1,M DO 30 J=1,M A(I,J)=A(I,J)+PP2(I,J)*WEI ADX(I,J)=ADX(I,J)+PP2(I,J)*DWEIX ADY(I,J)=ADY(I,J)+PP2(I,J)*DWEIY 30 CONTINUE C B(I,IN)=PB2(I)*WEI BDX(I,IN)=PB2(I)*DWEIX BDY(I,IN)=PB2(I)*DWEIY C 20 CONTINUE 10 CONTINUE C RETURN END C C C--------------------------------------------------------------------------- C BAS C C CALCULATE THE BASIS VECTOR FOR THE GIVEN NODE WHICH HAS THE INFLUENCE C AT THE GIVEN GAUSSIAN POINT. FIRST AND SECOND ORDER BASES ARE C PERMITTED. C--------------------------------------------------------------------------- C SUBROUTINE BAS(PB2,PP2,IN,ND,NN) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /GIPT/ X(2534),Y(2534) C DIMENSION ND(NN),PB2(6),PP2(6,6) C NI=ND(IN) C IF (M .EQ. 3) THEN PB2(1)=1.0 PB2(2)=X(NI) PB2(3)=Y(NI) ELSE IF (M .EQ. 6) THEN PB2(1)=1.0 PB2(2)=X(NI) PB2(3)=Y(NI) PB2(4)=X(NI)**2 PB2(5)=X(NI)*Y(NI) PB2(6)=Y(NI)**2 END IF C DO 10 I=1,M DO 20 J=1,M PP2(I,J)=PB2(I)*PB2(J) 20 CONTINUE 10 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C PVE C C CALCULATE THE BASIS VECTORS AND THEIR X,Y DERIVITIVES AT THE C GIVEN GAUSS POINT (XQ,YQ). FIRST AND SECOND ORDER BASES ARE C PERMITTED C------------------------------------------------------------------------ C SUBROUTINE PVE(PT,DPTX,DPTY,XQ,YQ) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M C DIMENSION PT(6),DPTX(6),DPTY(6) C IF (M .EQ. 3) THEN C PT(1)=1.0 PT(2)=XQ PT(3)=YQ C DPTX(1)=0.0 DPTX(2)=1.0 DPTX(3)=0.0 C DPTY(1)=0.0 DPTY(2)=0.0 DPTY(3)=1.0 C ELSE IF (M .EQ. 6) THEN C PT(1)=1.0 PT(2)=XQ PT(3)=YQ PT(4)=XQ**2 PT(5)=XQ*YQ PT(6)=YQ**2 C DPTX(1)=0.0 DPTX(2)=1.0 DPTX(3)=0.0 DPTX(4)=2*XQ DPTX(5)=YQ DPTX(6)=0.0 C DPTY(1)=0.0 DPTY(2)=0.0 DPTY(3)=1.0 DPTY(4)=0.0 DPTY(5)=XQ DPTY(6)=2*YQ C END IF C RETURN END C C C------------------------------------------------------------------------------- C SHAPE1 C C CALCULATE THE SHAPE FUNCTION (PHI) AND ITS X,Y DERIVATIVES C C WHERE C C PHI=P(X)A-1(X)B(X) C D(PHI)/DX=(DP/DX)(A-1)(B)+(P)[(A-1)(DB/DX)-(A-1)(DA/DX)(A-1)(B)] C D(PHI)/DY=(DP/DY)(A-1)(B)+(P)[(A-1)(DB/DY)-(A-1)(DA/DY)(A-1)(B)] C------------------------------------------------------------------------------- C SUBROUTINE SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M C DIMENSION IPIV(M) DIMENSION PHI(NN),DPHIX(NN),DPHIY(NN) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION A(M,M),B(M,NN),ADX(M,M),ADY(M,M),BDX(M,NN),BDY(M,NN) DIMENSION AINV(M,M),DXAINV(M,M),DYAINV(M,M),DAINVX(M,M), & DAINVY(M,M),AA(M,NN),AABDX(M,NN),AABDY(M,NN), & DAABX(M,NN),DAABY(M,NN),DPH1X(NN),DPH2X(NN),DPH1Y(NN), & DPH2Y(NN) C C CALL LAPACK SUBROUTINES C C AINV=A-1 C DO 15 I=1,M DO 16 J=1,M AINV(I,J)=0.0 16 CONTINUE AINV(I,I)=1.0 15 CONTINUE C CALL DGESV( M, M, A, M, IPIV, AINV, M, INFO ) C C DXAINV=(DA/DX)(A-1), DYAINV=(DA/DY)(A-1) C DO 10 I=1,M DO 20 J=1,M DXAINV(I,J)=0.0 DYAINV(I,J)=0.0 DO 30 K=1,M DXAINV(I,J)=DXAINV(I,J)+ADX(I,K)*AINV(K,J) DYAINV(I,J)=DYAINV(I,J)+ADY(I,K)*AINV(K,J) 30 CONTINUE 20 CONTINUE 10 CONTINUE C C DAINVX=-(A-1)(DA/DX)(A-1), DAINVY=-(A-1)(DA/DY)(A-1) C DO 40 I=1,M DO 50 J=1,M DAINVX(I,J)=0.0 DAINVY(I,J)=0.0 DO 60 K=1,M DAINVX(I,J)=DAINVX(I,J)-AINV(I,K)*DXAINV(K,J) DAINVY(I,J)=DAINVY(I,J)-AINV(I,K)*DYAINV(K,J) 60 CONTINUE 50 CONTINUE 40 CONTINUE C C AA=(A-1)(B) C AABDX=(A-1)(DB/DX), AABDY=(A-1)(DB/DY) C DO 70 I=1,M DO 80 J=1,NN AA(I,J)=0.0 AABDX(I,J)=0.0 AABDY(I,J)=0.0 DO 90 K=1,M AA(I,J)=AA(I,J)+AINV(I,K)*B(K,J) AABDX(I,J)=AABDX(I,J)+AINV(I,K)*BDX(K,J) AABDY(I,J)=AABDY(I,J)+AINV(I,K)*BDY(K,J) 90 CONTINUE 80 CONTINUE 70 CONTINUE C C DAABX=-(A-1)(DA/DX)(A-1)(B),DAABY=-(A-1)(DA/DY)(A-1)(B) C DO 100 I=1,M DO 110 J=1,NN DAABX(I,J)=0.0 DAABY(I,J)=0.0 DO 120 K=1,M DAABX(I,J)=DAABX(I,J)+DAINVX(I,K)*B(K,J) DAABY(I,J)=DAABY(I,J)+DAINVY(I,K)*B(K,J) 120 CONTINUE 110 CONTINUE 100 CONTINUE C C PHI=(P)(A-1)(B) C DPH1X=(DP/DX)(A-1)(B),DPH1Y=(DP/DY)(A-1)(B) C DPH2X=(A-1)(DB/DX)-(A-1)(DA/DX)(A-1)(B),DPH2Y=(A-1)(DB/DY)-(A-1)(DA/DY)(A-1)(B) C DPHIX=D(PHI)/DX=DPH1X+DPH2X C DPHIY=D(PHI)/DY=DPH1Y+DPH2Y C DO 130 I=1,NN PHI(I)=0.0 DPH1X(I)=0.0 DPH2X(I)=0.0 DPH1Y(I)=0.0 DPH2Y(I)=0.0 DO 140 J=1,M PHI(I)=PHI(I)+PT(J)*AA(J,I) DPH1X(I)=DPH1X(I)+DPTX(J)*AA(J,I) DPH2X(I)=DPH2X(I)+PT(J)*(AABDX(J,I)+DAABX(J,I)) DPH1Y(I)=DPH1Y(I)+DPTY(J)*AA(J,I) DPH2Y(I)=DPH2Y(I)+PT(J)*(AABDY(J,I)+DAABY(J,I)) 140 CONTINUE DPHIX(I)=DPH1X(I)+DPH2X(I) DPHIY(I)=DPH1Y(I)+DPH2Y(I) 130 CONTINUE C RETURN END C C C------------------------------------------------------------------------- C NSTIFF C C PERFORMS THE NUMERICAL INTEGRATION FOR THE STIFFNESS (SK) AND THE C MASS (SM) OVER THE NODAL DOMAIN C C WHERE C C SK=INTEGRAL{DVTI*D*BM} C SM=INTEGRAL(CT*PHI*VTI) C------------------------------------------------------------------------- C SUBROUTINE NSTIFF(SK,SM,IEL,NN,ND,VTI,STREST,PHI,DPHIX,DPHIY,WTX & ,WTY,DGJ) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /CON/ D(3,3) COMMON /NOTN/ NNODE COMMON /DEN/ CT C DIMENSION SK(2*NNODE,2*NNODE),SM(2*NNODE) DIMENSION ND(NN) DIMENSION SKK(2,2),BM(3,2),DBM(3,2),STREST(2,3) DIMENSION PHI(NN),DPHIX(NN),DPHIY(NN) C IKEL=2*IEL C DO 10 JN=1,NN C JND=2*ND(JN) C CALL BMATRX(BM,DPHIX,DPHIY,JN,NN) C DO 20 I=1,3 DO 30 J=1,2 DBM(I,J)=0.0 DO 40 K=1,3 DBM(I,J)=DBM(I,J)+D(I,K)*BM(K,J) 40 CONTINUE 30 CONTINUE 20 CONTINUE C DO 50 I=1,2 DO 60 J=1,2 SKK(I,J)=0.0 DO 70 K=1,3 SKK(I,J)=SKK(I,J)+STREST(I,K)*DBM(K,J) 70 CONTINUE 60 CONTINUE 50 CONTINUE C SM(IKEL-1)=SM(IKEL-1)+WTX*WTY*CT*PHI(JN)*VTI*DGJ SM(IKEL)=SM(IKEL-1) C SK(IKEL-1,JND-1)=SK(IKEL-1,JND-1)+WTX*WTY*SKK(1,1)*DGJ SK(IKEL-1,JND)=SK(IKEL-1,JND)+WTX*WTY*SKK(1,2)*DGJ SK(IKEL,JND-1)=SK(IKEL,JND-1)+WTX*WTY*SKK(2,1)*DGJ SK(IKEL,JND)=SK(IKEL,JND)+WTX*WTY*SKK(2,2)*DGJ C 10 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C BMATRX C C CALCULATE THE SHAPE FUNCTION DERIVATIVE MATRIX C------------------------------------------------------------------------ C SUBROUTINE BMATRX(BM,DPHIX,DPHIY,IN,NN) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION BM(3,2) DIMENSION DPHIX(NN),DPHIY(NN) C BM(1,1)=DPHIX(IN) BM(1,2)=0.0 BM(2,1)=0.0 BM(2,2)=DPHIY(IN) BM(3,1)=DPHIY(IN) BM(3,2)=DPHIX(IN) C RETURN END C C C----------------------------------------------------------------------- C BCAPLY C C PERFORM THE NUMERICAL INTEGRATION OF THE LOADING VECTOR (F) C AND THE BOUNDARY PORTION OF THE STIFFNESS MATRIX (SK) C C F=INTEGRAL{(VTI)(T_PRESCRIBED)+PENA(VTI)(NSS)(U_PRESCRIBED)} C C WHERE, C C PENA = PENALTY PARAMETER C VTI = TEST FUNCTION FOR NODE I AT THE GAUSS POINT C NSS = SPECIFIED DISPLACEMENT MATRIX C T_PRESCRIBED=PRESCRIBED TRACTIONS C U_PRESCRIBED=PRESCRIBED DISPLACEMENTS C----------------------------------------------------------------------- C SUBROUTINE BCAPLY(SK,F,POINT1,WEIGHT1) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /ROCCR/ RO(2534),CC(2534),R(2534) COMMON /GIPT/ X(2534),Y(2534) COMMON /NIPT/ SDSPX(2534),SDSPY(2534),SLDX(2534),SLDY(2534) COMMON /BOUND/ NBOUND,IBOUND(2534),INDBND(2534) COMMON /SPEC/ NDSPX(2534),NDSPY(2534),NLDX(2534),NLDY(2534) COMMON /GAU1/ NGL,NGLX,NGLY COMMON /PEN/ PENA C DIMENSION SK(2*NNODE,2*NNODE),F(2*NNODE) DIMENSION POINT1(9),WEIGHT1(9) DIMENSION ENDI(2,3),NSS(2,2) C LOGICAL SPECDX2,SPECDY2,SPECDX3,SPECDY3, & SPECLX2,SPECLY2,SPECLX3,SPECLY3 C C LOOPS OVER THE BOUNDARY NODES C DO 10 I=1,NBOUND C IEL=INDBND(I) IKEL=2*IEL C C FORWARD BOUDARY (I.E., BETWEEN I AND I+1 BOUNDARY NODES) C SPECDX3=NDSPX(IEL).LT.0.AND.NDSPX(IEL).NE.-3 SPECDY3=NDSPY(IEL).LT.0.AND.NDSPY(IEL).NE.-3 SPECLX3=NLDX(IEL).LT.0.AND.NLDX(IEL).NE.-3 SPECLY3=NLDY(IEL).LT.0.AND.NLDY(IEL).NE.-3 IF (I .EQ. NBOUND) THEN R1X = X(IEL) R1Y = Y(IEL) R2X = X(INDBND(1)) R2Y = Y(INDBND(1)) ELSE R1X = X(IEL) R1Y = Y(IEL) R2X = X(INDBND(I+1)) R2Y = Y(INDBND(I+1)) END IF C RMAG=DSQRT((R2X-R1X)**2+(R2Y-R1Y)**2) C C UNIT VECTOR IN DIRECTION OF THE BOUNDARY (ASSUMES STRAIGHT LINE C PATH C EX=(R2X-R1X)/RMAG EY=(R2Y-R1Y)/RMAG C C LOOP OVER THE GAUSS POINTS C DO 20 J=1,NGL C SKI=POINT1(J) WT=WEIGHT1(J) DGJ1=RO(IEL)/2.0 C C GLOBAL GAUSS POINT COOORDINATES C XB=R1X+0.5*RO(IEL)*(1+SKI)*EX YB=R1Y+0.5*RO(IEL)*(1+SKI)*EY C C GTHETA=DATAN2(YB,XB) C CXB=DCOS(GTHETA) C SYB=DSIN(GTHETA) C C UNIT OUTWARD NORMAL VECTOR C ENX=EY ENY=-EX C CALL NORMAL(ENDI,ENX,ENY) C C SPECIFIED DISPLACEMENT ENFORCEMENT MATRIX (FOR SPECIFIED X-DISP, C NSS(1,1)=1, FOR SPECIFIED Y-DISP NSS(2,2)=1 C NSS(1,1)=0 NSS(1,2)=0 NSS(2,1)=0 NSS(2,2)=0 C IF (SPECDX3) THEN UXSPD=SDSPX(IEL) NSS(1,1)=1 END IF C IF (SPECDY3) THEN UYSPD=SDSPY(IEL) NSS(2,2)=1 END IF C C DETERMINE LOADS TO ENFORCE THE ESSENTIAL BCS (PENALTY FUNCTION) C CALL GKCALC(SK,IEL,NSS,XB,YB,WT,VTI,DGJ1,ENDI) C IF (SPECDX3) THEN F(IKEL-1)=F(IKEL-1)+WT*PENA*VTI*UXSPD*DGJ1 END IF C IF (SPECDY3) THEN F(IKEL)=F(IKEL)+WT*PENA*VTI*UYSPD*DGJ1 END IF C C CALCULATE THE LOAD VECTOR FOR THE APPLIED LOADS (F=INTEGRAL{VTI*T}) C WHERE, VTI=TEST FUNCTION FOR NODE I AT THE GAUSS POINT C IF (SPECLX3) THEN F(IKEL-1)=F(IKEL-1)+WT*VTI*SLDX(IEL)*DGJ1 END IF C IF (SPECLY3) THEN F(IKEL)=F(IKEL)+WT*VTI*SLDY(IEL)*DGJ1 END IF C 20 CONTINUE C C BACKWARD BOUDARY (I.E., BETWEEN I-1 AND I BOUNDARY NODES) C SPECDX2=NDSPX(IEL).LT.0.AND.NDSPX(IEL).NE.-2 SPECDY2=NDSPY(IEL).LT.0.AND.NDSPY(IEL).NE.-2 SPECLX2=NLDX(IEL).LT.0.AND.NLDX(IEL).NE.-2 SPECLY2=NLDY(IEL).LT.0.AND.NLDY(IEL).NE.-2 C IF (I .EQ. 1) THEN R1X = X(IEL) R1Y = Y(IEL) R2X = X(INDBND(NBOUND)) R2Y = Y(INDBND(NBOUND)) ELSE R1X = X(IEL) R1Y = Y(IEL) R2X = X(INDBND(I-1)) R2Y = Y(INDBND(I-1)) END IF C RMAG=DSQRT((R2X-R1X)**2+(R2Y-R1Y)**2) C C UNIT VECTOR IN DIRECTION OF THE BOUNDARY (ASSUMES STRAIGHT LINE C PATH C EX=(R2X-R1X)/RMAG EY=(R2Y-R1Y)/RMAG C DO 30 J=1,NGL C SKI=POINT1(J) WT=WEIGHT1(J) DGJ1=RO(IEL)/2.0 C C GLOBAL GAUSS POINT COOORDINATES C XB=R1X+0.5*RO(IEL)*(1+SKI)*EX YB=R1Y+0.5*RO(IEL)*(1+SKI)*EY C C GTHETA=DATAN2(YB,XB) C CXB=DCOS(GTHETA) C SYB=DSIN(GTHETA) C C UNIT OUTWARD NORMAL VECTOR C ENX=-EY ENY=EX C CALL NORMAL(ENDI,ENX,ENY) C C SPECIFIED DISPLACEMENT ENFORCEMENT MATRIX (FOR SPECIFIED X-DISP, C NSS(1,1)=1, FOR SPECIFIED Y-DISP NSS(2,2)=1 C NSS(1,1)=0 NSS(1,2)=0 NSS(2,1)=0 NSS(2,2)=0 C IF (SPECDX2) THEN UXSPD=SDSPX(IEL) NSS(1,1)=1 END IF C IF (SPECDY2) THEN UYSPD=SDSPY(IEL) NSS(2,2)=1 END IF C C DETERMINE LOADS TO ENFORCE THE ESSENTIAL BCS (PENALTY FUNCTION) C CALL GKCALC(SK,IEL,NSS,XB,YB,WT,VTI,DGJ1,ENDI) C IF (SPECDX2) THEN F(IKEL-1)=F(IKEL-1)+WT*PENA*VTI*UXSPD*DGJ1 END IF C IF (SPECDY2) THEN F(IKEL)=F(IKEL)+WT*PENA*VTI*UYSPD*DGJ1 END IF C C CALCULATE THE LOAD VECTOR FOR THE APPLIED LOADS (F=INTEGRAL{VTI*T}) C WHERE, VTI=TEST FUNCTION FOR NODE I AT THE GAUSS POINT C IF (SPECLX2) THEN F(IKEL-1)=F(IKEL-1)+WT*VTI*SLDX(IEL)*DGJ1 END IF C IF (SPECLY2) THEN F(IKEL)=F(IKEL)+WT*VTI*SLDY(IEL)*DGJ1 END IF C 30 CONTINUE 10 CONTINUE C RETURN END C C C----------------------------------------------------------------------- C GKCALC C C PERFORM THE NUMERICAL INTEGRATION OF THE BOUNDARY PORTION OF THE C STIFFNESS MATRIX (SK) C C SK=INTEGRAL{(PENA)(NSS)(VTI)(PHI)-(VTI)(ENDI)(D)(BM)(NSS)} C C WHERE, C C PENA = PENALTY PARAMETER (STIFFNESS) C VTI = TEST FUNCTION FOR NODE I AT THE GAUSS POINT C PHI = SHAPE FUNCTION C ENDI = UNIT OUTWARD NORMAL MATRIX C D = CONSTITUTIVE MATRIX C NSS = SPECIFIED DISPLACEMENT MATRIX C BM = SHAPE FUNCTION DERIVITIVE MATRIX C----------------------------------------------------------------------- C SUBROUTINE GKCALC(SK,IEL,NSS,XB,YB,WT,VTI,DGJ1,ENDI) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /CON/ D(3,3) COMMON /PEN/ PENA COMMON /ROCCR/ RO(2534),CC(2534),R(2534) C DIMENSION SK(2*NNODE,2*NNODE) DIMENSION BM(3,2),ENDI(2,3),STREST(2,3) DIMENSION NSS(2,2),SKU1(2,2),SKU2(2,2),VND(2,3),VNDB(2,2) DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) C C CALCULATE ALL THE REQUIRED VARIABLES C IKEL=2*IEL C CALL WCOUNT(XB,YB,NN,ND,W,DWX,DWY) C CALL TESFUN(XB,YB,IEL,VTI,STREST) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XB,YB) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C C INTEGRATE OVER ALL OF THE NODES IN THE DOMAIN OF DEFINITION FOR C THE GIVEN GAUSS POINT C DO 20 IN=1,NN C IND=2*ND(IN) C CALL BMATRX(BM,DPHIX,DPHIY,IN,NN) C C SKU1=(NSS)(VTI)(PHI) C DO 30 I=1,2 DO 40 J=1,2 SKU1(I,J)=NSS(I,J)*VTI*PHI(IN) 40 CONTINUE 30 CONTINUE C C SKU2=(VTI)(ENDI)(D)(BM)(NSS) C DO 50 I=1,2 DO 60 J=1,3 VND(I,J)=0.0 DO 70 K=1,3 VND(I,J)=VND(I,J)+VTI*ENDI(I,K)*D(K,J) 70 CONTINUE 60 CONTINUE 50 CONTINUE C DO 80 I=1,2 DO 90 J=1,2 VNDB(I,J)=0.0 DO 100 K=1,3 VNDB(I,J)=VNDB(I,J)+VND(I,K)*BM(K,J) 100 CONTINUE 90 CONTINUE 80 CONTINUE C DO 110 I=1,2 DO 120 J=1,2 SKU2(I,J)=0.0 DO 130 K=1,2 SKU2(I,J)=SKU2(I,J)+VNDB(I,K)*NSS(K,J) 130 CONTINUE 120 CONTINUE 110 CONTINUE C SK(IKEL-1,IND-1)=SK(IKEL-1,IND-1)+WT*PENA*SKU1(1,1)*DGJ1 & -WT*SKU2(1,1)*DGJ1 SK(IKEL-1,IND)=SK(IKEL-1,IND)+WT*PENA*SKU1(1,2)*DGJ1 & -WT*SKU2(1,2)*DGJ1 SK(IKEL,IND-1)=SK(IKEL,IND-1)+WT*PENA*SKU1(2,1)*DGJ1 & -WT*SKU2(2,1)*DGJ1 SK(IKEL,IND)=SK(IKEL,IND)+WT*PENA*SKU1(2,2)*DGJ1 & -WT*SKU2(2,2)*DGJ1 C 20 CONTINUE C RETURN END C C C------------------------------------------------------------------------ C NORMAL C C CREATE THE UNIT OUTWARD NORMAL MATRIX (ENDI) C------------------------------------------------------------------------ C SUBROUTINE NORMAL(ENDI,ENX,ENY) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION ENDI(2,3) C ENDI(1,1)=ENX ENDI(1,2)=0.0 ENDI(1,3)=ENY ENDI(2,1)=0.0 ENDI(2,2)=ENY ENDI(2,3)=ENX C RETURN END C C C------------------------------------------------------------------------ C IACCE C C SOLVE FOR THE INITIAL ACCELERATION C------------------------------------------------------------------------ C SUBROUTINE IACCE(SM,F) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION SM(2*NNODE),F(2*NNODE) C NEQ=2*NNODE C DO 10 I=1,NEQ GLA(I)=F(I)/SM(I) 10 CONTINUE C RETURN END C C C--------------------------------------------------------------------- C TIMESKM (NEWMARK FAMILY OF METHODS) C C DETERMINE THE EFFECTIVE STIFFNESS MATRIX (SKM) C C SKM=SK+A3*SM C--------------------------------------------------------------------- C SUBROUTINE TIMESKM(SKM,SK,SM) C IMPLICIT REAL*8(A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /TS/ DT,A1,A2,A3,A4,A5 C DIMENSION SKM(2*NNODE,2*NNODE),SK(2*NNODE,2*NNODE),SM(2*NNODE) C C THE NEWMARK INTERGRATION SCHEME FOR HYPERBOLIC EQUATIONS C NEQ=2*NNODE C DO 10 I=1,NEQ SKM(I,I)=A3*SM(I) DO 20 J=1,NEQ SKM(I,J)=SKM(I,J)+SK(I,J) 20 CONTINUE 10 CONTINUE C RETURN END C C C--------------------------------------------------------------------- C INVEKM (NEWMARK FAMILY OF METHODS) C C SOLVE FOR THE INVERSE OF THE EFFECTIVE STIFFNESS MATRIX (SKM) C C SK=INVERSE(SKM) C--------------------------------------------------------------------- C SUBROUTINE INVEKM(SKM,SK) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE C DIMENSION IPIV(2*NNODE) DIMENSION SKM(2*NNODE,2*NNODE),SK(2*NNODE,2*NNODE) C NEQ=2*NNODE C C CALL LAPACK SUBROUTINES C DO 10 I=1,NEQ DO 20 J=1,NEQ SK(I,J)=0.0 20 CONTINUE SK(I,I)=1.0 10 CONTINUE C CALL DGESV( NEQ, NEQ, SKM, NEQ, IPIV, SK, NEQ, INFO ) C RETURN END C C C---------------------------------------------------------------------- C TIMEFM (NEWMARK FAMILY OF METHODS) C C DETERMINE THE EFFECTIVE LOADING VECTOR (FM) C C FM=F+SM*(A3*GLU+A4*GLV+A5*GLA) C---------------------------------------------------------------------- C SUBROUTINE TIMEFM(FM,F,SM) C IMPLICIT REAL*8(A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /TS/ DT,A1,A2,A3,A4,A5 COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION FM(2*NNODE),F(2*NNODE),SM(2*NNODE) C C THE NEWMARK INTERGRATION SCHEME FOR HYPERBOLIC EQUATIONS C NEQ=2*NNODE C DO 10 I=1,NEQ C FM(I)=F(I)+SM(I)*(A3*GLU(I)+A4*GLV(I)+A5*GLA(I)) C 10 CONTINUE C RETURN END C C C---------------------------------------------------------------------- C SOLVE (NEWMARK FAMILY OF METHODS) C C SOLVE FOR THE DISPLACEMENTS AT THE (N+1) TIME STEP C---------------------------------------------------------------------- C SUBROUTINE SOLVE(SK,FM,FSOL) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE C DIMENSION SK(2*NNODE,2*NNODE),FM(2*NNODE),FSOL(2*NNODE) C NEQ=2*NNODE C DO 10 I=1,NEQ DO 20 J=1,NEQ FSOL(I)=FSOL(I)+SK(I,J)*FM(J) 20 CONTINUE 10 CONTINUE C RETURN END C C C--------------------------------------------------------------------------- C NVEAC (NEWMARK FAMILY OF METHODS) C C OBTAIN THE VELOCITY AND ACCELERATION AT TIME STEP (N+1) C--------------------------------------------------------------------------- C SUBROUTINE NVEAC(FSOL) C IMPLICIT REAL*8(A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /TS/ DT,A1,A2,A3,A4,A5 COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION FSOL(2*NNODE) C NEQ=2*NNODE C DO 10 I=1,NEQ GLU(I)=A3*(FSOL(I)-GLU(I))-A4*GLV(I)-A5*GLA(I) GLV(I)=GLV(I)+A1*GLU(I)+A2*GLA(I) GLA(I)=GLU(I) GLU(I)=FSOL(I) 10 CONTINUE C RETURN END C C C---------------------------------------------------------------------------- C TIMER (THE CENTRAL DIFFERENCE METHOD) C C OBTAIN THE DISPLACEMENT, VELOCITY, AND ACCELERATION AT TIME STEP (N+1) C---------------------------------------------------------------------------- C SUBROUTINE TIMER(FM,FSOL,SK,SM,F) C IMPLICIT REAL*8(A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /TS/ DT,A1,A2,A3,A4,A5 COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION FM(2*NNODE),FSOL(2*NNODE) DIMENSION SK(2*NNODE,2*NNODE),SM(2*NNODE),F(2*NNODE) C NEQ=2*NNODE C DO 10 I=1,NEQ GLU(I)=GLU(I)+DT*GLV(I)+(DT**2)/2*GLA(I) 10 CONTINUE C DO 20 I=1,NEQ C SUM=0.0 C DO 30 J=1,NEQ SUM=SUM+SK(I,J)*GLU(J) 30 CONTINUE C FM(I)=F(I)-SUM C FSOL(I)=FM(I)/SM(I) C GLV(I)=GLV(I)+DT/2*(FSOL(I)+GLA(I)) C GLA(I)=FSOL(I) C 20 CONTINUE C RETURN END C C C-------------------------------------------------------------------------------- C OUTP C C FOR THE OUTPUT OF THE TIME HISTORY OF DISPLACEMENTS AND STRESSES AT C GIVEN NODES C C NRSE,NSG(I) : NUMBER OF THE NODES DIRECTLY AHEAD OF THE CRACK TIP AND THEIR C NODE NO. (RESPECTIVELY) C-------------------------------------------------------------------------------- C SUBROUTINE OUTP(ICT,STIST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /NOUP/ NRSE,NSG(6) C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) C DO 10 ICOUNT=1,NRSE C ISD=NSG(ICOUNT) C XN=X(ISD) YN=Y(ISD) C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIGT1(ICT,STIST,ICOUNT,XN,YN,NN,ND,PHI,DPHIX,DPHIY) C 10 CONTINUE C RETURN END C C C----------------------------------------------------------------------------- C DISIGT1 C C PRINT THE DISPLACEMENTS AND STRESSES C----------------------------------------------------------------------------- C SUBROUTINE DISIGT1(ICT,STIST,ICOUNT,XN,YN,NN,ND,PHI,DPHIX,DPHIY) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /CON/ D(3,3) COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION PHI(NN),DPHIX(NN),DPHIY(NN),ND(NN) DIMENSION STRAIN(3),SIGMA(3) C C OPEN (UNIT=21,FILE='disp1.dat',STATUS='UNKNOWN') C OPEN (UNIT=22,FILE='disp2.dat',STATUS='UNKNOWN') C OPEN (UNIT=23,FILE='disp3.dat',STATUS='UNKNOWN') C OPEN (UNIT=24,FILE='disp4.dat',STATUS='UNKNOWN') C OPEN (UNIT=25,FILE='disp5.dat',STATUS='UNKNOWN') C OPEN (UNIT=26,FILE='disp6.dat',STATUS='UNKNOWN') OPEN (UNIT=31,FILE='stress1.dat',STATUS='UNKNOWN') OPEN (UNIT=32,FILE='stress2.dat',STATUS='UNKNOWN') OPEN (UNIT=33,FILE='stress3.dat',STATUS='UNKNOWN') C OPEN (UNIT=34,FILE='stress4.dat',STATUS='UNKNOWN') C OPEN (UNIT=35,FILE='stress5.dat',STATUS='UNKNOWN') C OPEN (UNIT=36,FILE='stress6.dat',STATUS='UNKNOWN') C UI=0.0 VI=0.0 DUIX=0.0 DVIX=0.0 DVIY=0.0 DUIY=0.0 C DO 10 IN=1,NN C IND=2*ND(IN) C UI=UI+PHI(IN)*GLU(IND-1) VI=VI+PHI(IN)*GLU(IND) DUIX=DUIX+DPHIX(IN)*GLU(IND-1) DVIX=DVIX+DPHIX(IN)*GLU(IND) DVIY=DVIY+DPHIY(IN)*GLU(IND) DUIY=DUIY+DPHIY(IN)*GLU(IND-1) C 10 CONTINUE C STRAIN(1)=DUIX STRAIN(2)=DVIY STRAIN(3)=DVIX+DUIY C DO 20 I=1,3 SIGMA(I)=0.0 DO 30 J=1,3 SIGMA(I)=SIGMA(I)+D(I,J)*STRAIN(J) 30 CONTINUE 20 CONTINUE C C LOUTD=20 LOUTS=30 C C LOUTD=LOUTD+ICOUNT LOUTS=LOUTS+ICOUNT C IF (ICT .EQ. 1) THEN C C WRITE(LOUTD,'(A)') ' TIME X-COORD Y-COORD C & X-DISP Y-DISP' C WRITE(LOUTD,100) C WRITE(LOUTS,'(A)') ' TIME X-COORD Y-COORD & STRESS-X STRESS-Y STRESS-XY' WRITE(LOUTS,100) C END IF C C PRINT OUT THE DISPLACEMENTS AND STRESSES C C WRITE(LOUTD,110) STIST,XN,YN,UI,VI C WRITE(LOUTS,120) STIST,XN,YN,(SIGMA(I),I=1,3) C 100 FORMAT(1H ) C110 FORMAT(1X,3E15.6,2E16.6) 120 FORMAT(1X,3E15.6,3E16.6) C RETURN END C C C-------------------------------------------------------------------------------- C OUTALL C C FOR THE OUTPUT OF DISPLACEMENTS AND STRESSES OVER THE WHOLE DOMAIN AT THE C SPECIFIED TIME STEP (NOUPT) C-------------------------------------------------------------------------------- C SUBROUTINE OUTALL(ICT,STIST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /LOAD/ XLOAD C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) DIMENSION UU(NNODE),VV(NNODE),STRX(NNODE),STRY(NNODE),STRXY(NNODE) C OPEN (UNIT=46,FILE='dispall.dat',STATUS='UNKNOWN') OPEN (UNIT=47,FILE='strall.dat',STATUS='UNKNOWN') OPEN (UNIT=48,FILE='prins.dat',STATUS='UNKNOWN') C DO 10 ICOUNT=1,NNODE C XN=X(ICOUNT) YN=Y(ICOUNT) C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND,PHI,DPHIX,DPHIY) C 10 CONTINUE C C FOR DISPLACEMENTS C WRITE(46,'(A)')' TIME STEP #, TIME' WRITE(46,*) ICT,STIST WRITE(46,100) C WRITE(46,'(A)')' NODE# X-COORD Y-COORD &X-DISP Y-DISP' WRITE(46,100) C DO 20 ID=1,NNODE C WRITE(46,110) ID,X(ID),Y(ID),UU(ID),VV(ID) C 20 CONTINUE C C FOR X, Y STRESS COMPONENTS C WRITE(47,'(A)')' TIME STEP #, TIME' WRITE(47,*) ICT,STIST WRITE(47,100) C WRITE(47,'(A)')' NODE# X-COORD Y-COORD &STRESS-X STRESS-Y STRESS-XY' WRITE(47,100) C C FOR PRINCIPAL STRESS COMPONENTS C WRITE(48,'(A)')' TIME STEP #, TIME' WRITE(48,*) ICT,STIST WRITE(48,100) C WRITE(48,'(A)')' NODE# X-COORD Y-COORD &MAX-PRINP MINS-PRINP MAX-SHEAR' WRITE(48,100) C DO 30 IS=1,NNODE C WRITE(47,120) IS,X(IS),Y(IS),STRX(IS),STRY(IS),STRXY(IS) C SSA=(STRX(IS)-STRY(IS))/2.0 C C PRINCIPAL STRESSES NORMALIZED BY THE APPLIED TRACTION C SMAX=((STRX(IS)+STRY(IS))/2.0+DSQRT(SSA**2+STRXY(IS)**2))/XLOAD SMIN=((STRX(IS)+STRY(IS))/2.0-DSQRT(SSA**2+STRXY(IS)**2))/XLOAD TAUM=(DSQRT(SSA**2+STRXY(IS)**2))/XLOAD C WRITE(48,120) IS,X(IS),Y(IS),SMAX,SMIN,TAUM C 30 CONTINUE C 100 FORMAT(1H ) 110 FORMAT(1X,I4,4E17.7) 120 FORMAT(1X,I4,5E16.6) C RETURN END C C C-------------------------------------------------------------------------------- C OUTPA C C FOR THE OUTPUT OF DISPLACEMENTS AND STRESSES ALONG THE LINE AHEAD OF THE C CRACK TIP AT THE SPECIFIED TIME STEP (NOUPT) C C SL, NINP : RADIAL DISTANCE DIRECTLY AHEAD OF THE CRACK TIP AND NUMBER OF C POINTS FOR OUTPUT (RESPECTIVELY) C-------------------------------------------------------------------------------- C SUBROUTINE OUTPA(ICT,STIST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /CRACK/ CX1,CY1,CX2,CY2 COMMON /OUPD/ SL,NINP C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) DIMENSION UU(NNODE),VV(NNODE),STRX(NNODE),STRY(NNODE),STRXY(NNODE) DIMENSION RADI(NINP) C C DATA PI /3.141592653/ C OPEN (UNIT=56,FILE='dispfar.dat',STATUS='UNKNOWN') OPEN (UNIT=57,FILE='strfar.dat',STATUS='UNKNOWN') C DX=SL/(NINP-1) C DO 10 ICOUNT=1,NINP C RADI(ICOUNT)=(ICOUNT-1)*DX XN=CX2+RADI(ICOUNT) YN=CY2 C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND,PHI,DPHIX,DPHIY) C 10 CONTINUE C C FOR DISPLACEMENTS C WRITE(56,'(A)')' TIME STEP #, TIME' WRITE(56,*) ICT,STIST WRITE(56,100) C WRITE(56,'(A)')' NUM# RADI-DISTANCE X-DISP & Y-DISP' WRITE(56,100) C DO 20 ID=1,NINP C WRITE(56,110) ID,RADI(ID),UU(ID),VV(ID) C 20 CONTINUE C C FOR STRESSES C WRITE(57,'(A)')' TIME STEP #, TIME' WRITE(57,*) ICT,STIST WRITE(57,100) C WRITE(57,'(A)')' NUM# RADI-DISTANCE STRESS-X STRESS-Y & STRESS-XY' WRITE(57,100) C DO 30 IS=1,NINP C WRITE(57,120) IS,RADI(IS),STRX(IS),STRY(IS),STRXY(IS) C 30 CONTINUE C 100 FORMAT(1H ) 110 FORMAT(1X,I4,3E18.8) 120 FORMAT(1X,I4,4E15.6) C RETURN END C C C-------------------------------------------------------------------------------- C DEFO C C FOR THE OUTPUT OF THE DEFORMED SHAPE OF THE DOMAIN AT THE SPECIFIED TIME C STEP (NOUPT) C C NBOUND : TOTAL NUMBER OF BOUNDARY NODES C INDBND(I) : THE NODE NUMBER FOR THE Ith NODE ALONG THE BOUNDARY C IBOUND(I*) : THE Ith NODE CORRESPONDING TO THE NODE NUMBER I* C-------------------------------------------------------------------------------- C SUBROUTINE DEFO(ICT,STIST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /GIPT/ X(2534),Y(2534) COMMON /BOUND/ NBOUND,IBOUND(2534),INDBND(2534) C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) DIMENSION UU(NNODE),VV(NNODE),STRX(NNODE),STRY(NNODE),STRXY(NNODE) DIMENSION DEFOX(NBOUND),DEFOY(NBOUND) C OPEN (UNIT=59,FILE='deform3.dat',STATUS='UNKNOWN') C DO 10 ICOUNT=1,NBOUND C XN=X(INDBND(ICOUNT)) YN=Y(INDBND(ICOUNT)) C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND,PHI,DPHIX,DPHIY) C 10 CONTINUE C C FOR BOUNDARY DEFORMATION C WRITE(59,'(A)')' TIME STEP #, TIME' WRITE(59,*) ICT,STIST WRITE(59,100) C WRITE(59,'(A)')'NODE# X-COORD Y-COORD X-DISP &Y-DISP DEFOR-X DEFOR-Y' WRITE(59,100) C DO 20 ID=1,NBOUND C DEFOX(ID)=X(INDBND(ID))+UU(ID) DEFOY(ID)=Y(INDBND(ID))+VV(ID) C WRITE(59,110) INDBND(ID),X(INDBND(ID)),Y(INDBND(ID)) & ,UU(ID),VV(ID),DEFOX(ID),DEFOY(ID) C 20 CONTINUE C 100 FORMAT(1H ) 110 FORMAT(1X,I4,6E13.5) C RETURN END C C C----------------------------------------------------------------------------- C DISIG C C PRINT THE DISPLACEMENTS AND STRESSES C----------------------------------------------------------------------------- C SUBROUTINE DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND, & PHI,DPHIX,DPHIY) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /CON/ D(3,3) COMMON /GTF/ GLU(5068),GLV(5068),GLA(5068) C DIMENSION PHI(NN),DPHIX(NN),DPHIY(NN),ND(NN) DIMENSION STRAIN(3),SIGMA(3) DIMENSION UU(NNODE),VV(NNODE),STRX(NNODE),STRY(NNODE),STRXY(NNODE) C UI=0.0 VI=0.0 DUIX=0.0 DVIX=0.0 DVIY=0.0 DUIY=0.0 C DO 10 IN=1,NN C IND=2*ND(IN) C UI=UI+PHI(IN)*GLU(IND-1) VI=VI+PHI(IN)*GLU(IND) DUIX=DUIX+DPHIX(IN)*GLU(IND-1) DVIX=DVIX+DPHIX(IN)*GLU(IND) DVIY=DVIY+DPHIY(IN)*GLU(IND) DUIY=DUIY+DPHIY(IN)*GLU(IND-1) C 10 CONTINUE C UU(ICOUNT)=UI VV(ICOUNT)=VI C STRAIN(1)=DUIX STRAIN(2)=DVIY STRAIN(3)=DVIX+DUIY C DO 20 I=1,3 SIGMA(I)=0.0 DO 30 J=1,3 SIGMA(I)=SIGMA(I)+D(I,J)*STRAIN(J) 30 CONTINUE 20 CONTINUE C STRX(ICOUNT)=SIGMA(1) STRY(ICOUNT)=SIGMA(2) STRXY(ICOUNT)=SIGMA(3) C RETURN END C C C------------------------------------------------------------------------------- C KCOMP C C COMPUTE THE STRESS INTENSITY FACTORS WITH THE NEAR-TIP STRESS FIELDS BY C USING THE LEAST SQUARES METHOD C C NEX : NUMBER OF POINTS FOR COMPUTING STRESS INTENSITY FACTOR ALONG A RANGE C (R) NEAR THE CRACK TIP C XK1A,XK1B,XK2A,XK2B : XK1A,XK1B = THE STARTING AND ENDING POINTS FOR THE C FIRST R RANGE C XK2A,XK2B = THE STARTING AND ENDING POINTS FOR THE C SECOND R RANGE C XLOAD : THE TRACTION FORCES (ASSUMED TO BE UNIFORM) C------------------------------------------------------------------------------- C SUBROUTINE KCOMP(ICT,STIST) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MBO/ M COMMON /NOTN/ NNODE COMMON /OUTK/ NEX COMMON /RANGE/ XK1A,XK1B,XK2A,XK2B COMMON /LOAD/ XLOAD COMMON /CRACK/ CX1,CY1,CX2,CY2 COMMON /OUTT/ NOUPT C DIMENSION ND(NNODE),W(NNODE),DWX(NNODE),DWY(NNODE) DIMENSION A(M,M),B(M,NNODE),ADX(M,M),ADY(M,M),BDX(M,NNODE), & BDY(M,NNODE) DIMENSION PT(6),DPTX(6),DPTY(6) DIMENSION PHI(NNODE),DPHIX(NNODE),DPHIY(NNODE) DIMENSION UU(NNODE),VV(NNODE),STRX(NNODE),STRY(NNODE),STRXY(NNODE) DIMENSION RADX(NEX) C DATA PI /3.141592653/ C OPEN (UNIT=76,FILE='k1k2.dat',STATUS='UNKNOWN') OPEN (UNIT=77,FILE='slope.dat',STATUS='UNKNOWN') OPEN (UNIT=78,FILE='tstep1.dat',STATUS='UNKNOWN') C C TO COMPUTE K1 FORM SIGMAX C C FOR THE FIRST R RANGE C DX=(XK1B-XK1A)/(NEX-1) C DO 10 ICOUNT=1,NEX C RADX(ICOUNT)=(ICOUNT-1)*DX+(XK1A-CX2) XN=CX2+RADX(ICOUNT) YN=CY2 C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND,PHI,DPHIX,DPHIY) C 10 CONTINUE C C CALL THE SUBROUTINE FOR THE LEAST SQUARES METHOD C CALL REGRE1(SKX,SLOP1,RADX,STRX) C C FOR THE NORMALIZED K1 (K1=K1/XLOAD*SQRT(PI*CX2)) C SKX=SKX/(XLOAD*DSQRT(PI*CX2)) C C TO COMPUTE K1 FORM SIGMAY C C FOR THE SECOND R RANGE C DX=(XK2B-XK2A)/(NEX-1) C DO 20 ICOUNT=1,NEX C RADX(ICOUNT)=(ICOUNT-1)*DX+(XK2A-CX2) XN=CX2+RADX(ICOUNT) YN=CY2 C CALL WCOUNT(XN,YN,NN,ND,W,DWX,DWY) C CALL INIAB(A,B,ADX,ADY,BDX,BDY,NN) C CALL ABM(A,B,ADX,ADY,BDX,BDY,NN,ND,W,DWX,DWY) C CALL PVE(PT,DPTX,DPTY,XN,YN) C CALL SHAPE1(PHI,DPHIX,DPHIY,A,ADX,ADY,B,BDX,BDY,PT,DPTX,DPTY & ,NN) C CALL DISIG(UU,VV,STRX,STRY,STRXY,ICOUNT,NN,ND,PHI,DPHIX,DPHIY) C 20 CONTINUE C C CALL THE SUBROUTINE FOR THE LEAST SQUARES METHOD C CALL REGRE2(SKY,SLOP2,RADX,STRY) C C FOR THE NORMALIZED K1 (K1=K1/XLOAD*SQRT(PI*CX2)) C SKY=SKY/(XLOAD*DSQRT(PI*CX2)) C C FOR THE OUTPUT OF THE K1 FORM SIGMAX AND SIGMAY RESPECTIVELY C IF (ICT .EQ. 1) THEN C WRITE(76,'(A)')'TIME STEP# TIME KX &KY' WRITE(76,100) C C FOR THE OUTPUT OF THE ORDER OF SINGULARITY FORM SIGMAX AND SIGMAY RESPECTIVELY C WRITE(77,'(A)')'TIME STEP# TIME SLOPEX &SLOPEY' WRITE(77,100) C END IF C WRITE(76,110) ICT,STIST,SKX,SKY C WRITE(77,110) ICT,STIST,SLOP1,SLOP2 C C OUTPUT THE STRESSES AT THE TIME STEP (NOUPT) C IF (ICT .EQ. NOUPT) THEN C WRITE(78,'(A)')'TIME STEP #, TIME' WRITE(78,*) ICT,STIST WRITE(78,100) C WRITE(78,'(A)')' NUM# RADI-DISTANCE STRESS-X & STRESS-Y STRESS-XY' WRITE(78,100) C DO 30 IS=1,NEX C WRITE(78,120) IS,RADX(IS),STRX(IS),STRY(IS),STRXY(IS) C 30 CONTINUE C END IF C 100 FORMAT(1H ) 110 FORMAT(1X,I5,3E16.6) 120 FORMAT(1X,I4,4E15.6) C RETURN END C C C------------------------------------------------------------------------------ C REGRE1 C C FITTING A STRAIGHT LINE THROUGH THE DATA OF STRESSES VS. RADIAL DISTANCES C (RADX) ON THE LOG SCALE BY THE LEAST SQUARES METHOD C------------------------------------------------------------------------------ C SUBROUTINE REGRE1(SK1,SLOP1,RADX,STRX) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /OUTK/ NEX C DIMENSION RADX(NEX),RLX(NEX) DIMENSION STRX(NNODE),STK1(NEX) DIMENSION AX(2,2),BX(2) DIMENSION IPIV(2) C DATA PI /3.141592653/ C DO 10 I=1,2 BX(I)=0.0 DO 20 J=1,2 AX(I,J)=0.0 20 CONTINUE 10 CONTINUE C DO 30 I=1,NEX C RLX(I)=DLOG(DABS(RADX(I))) STK1(I)=DLOG(DABS(STRX(I))) C AX(1,1)=AX(1,1)+1.0 AX(1,2)=AX(1,2)+RLX(I) AX(2,1)=AX(1,2) AX(2,2)=AX(2,2)+RLX(I)**2 C BX(1)=BX(1)+STK1(I) BX(2)=BX(2)+RLX(I)*STK1(I) C 30 CONTINUE C NRHS=1 C CALL DGESV(2,NRHS,AX,2,IPIV,BX,2,INFO) C C BX(1) IS THE INTERCEPTION OF Y AXIS IN THE LOG SCALE C BX(2) IS THE ORDER OF SINGULARITY C SK1=DSQRT(2*PI)*DEXP(BX(1)) SLOP1=BX(2) C RETURN END C C C------------------------------------------------------------------------------ C REGRE2 C C FITTING A STRAIGHT LINE THROUGH THE DATA OF STRESSES VS. RADIAL DISTANCES C (RADX) ON THE LOG SCALE BY THE LEAST SQUARES METHOD C------------------------------------------------------------------------------ C SUBROUTINE REGRE2(SK2,SLOP2,RADX,STRY) C IMPLICIT REAL*8 (A-H,O-Z) C COMMON /NOTN/ NNODE COMMON /OUTK/ NEX C DIMENSION RADX(NEX),RLX(NEX) DIMENSION STRY(NNODE),STK2(NEX) DIMENSION AX(2,2),BX(2) DIMENSION IPIV(2) C DATA PI /3.141592653/ C DO 10 I=1,2 BX(I)=0.0 DO 20 J=1,2 AX(I,J)=0.0 20 CONTINUE 10 CONTINUE C DO 30 I=1,NEX C RLX(I)=DLOG(DABS(RADX(I))) STK2(I)=DLOG(DABS(STRY(I))) C AX(1,1)=AX(1,1)+1.0 AX(1,2)=AX(1,2)+RLX(I) AX(2,1)=AX(1,2) AX(2,2)=AX(2,2)+RLX(I)**2 C BX(1)=BX(1)+STK2(I) BX(2)=BX(2)+RLX(I)*STK2(I) C 30 CONTINUE C NRHS=1 C CALL DGESV(2,NRHS,AX,2,IPIV,BX,2,INFO) C C BX(1) IS THE INTERCEPTION OF Y AXIS IN THE LOG SCALE C BX(2) IS THE ORDER OF SINGULARITY C SK2=DSQRT(2*PI)*DEXP(BX(1)) SLOP2=BX(2) C RETURN END C C * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module SORT2 from package NAPACK. * Retrieved from NETLIB on Tue Mar 28 01:36:29 2000. * ====================================================================== C C ________________________________________________________ C | | C | SORT AN ARRAY IN INCREASING ORDER | C | | C | INPUT: | C | | C | X --ARRAY OF NUMBERS | C | | C | W --WORKING ARRAY (LENGTH AT LEAST N) | C | | C | N --NUMBER OF ARRAY ELEMENTS TO SORT | C | | C | OUTPUT: | C | | C | X --ORIGINAL ARRAY | C | | C | Y --INDICES OF X GIVING INCREASING ORDER | C |________________________________________________________| C SUBROUTINE SORT2(X,Y,DW,N) C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 Y(1) DIMENSION DW(1),X(1) INTEGER I,J,K,L,M,N,P,Q C I = 1 10 K = I 20 J = I Y(I) = I I = I + 1 IF ( J .EQ. N ) GOTO 30 IF ( X(I) .GE. X(J) ) GOTO 20 DW(K) = I GOTO 10 30 IF ( K .EQ. 1 ) RETURN DW(K) = N + 1 40 M = 1 L = 1 50 I = L IF ( I .GT. N ) GOTO 120 P = Y(I) S = X(P) J = DW(I) K = J IF ( J .GT. N ) GOTO 100 Q = Y(J) T = X(Q) L = DW(J) Y(I) = L 60 IF ( S .GT. T ) GOTO 70 DW(M) = P M = M + 1 I = I + 1 IF ( I .EQ. K ) GOTO 80 P = Y(I) S = X(P) GOTO 60 70 DW(M)= Q M = M + 1 J = J + 1 IF ( J .EQ. L ) GOTO 110 Q = Y(J) T = X(Q) GOTO 60 80 DW(M) = Q K = M + L - J I = J - M 90 M = M + 1 IF ( M .EQ. K ) GOTO 50 DW(M) = Y(M+I) GOTO 90 100 Y(I) = J L = J 110 DW(M) = P K = M + K - I I = I - M GOTO 90 120 I = 1 130 K = I J = Y(I) 140 Y(I) = DW(I) I = I + 1 IF ( I .LT. J ) GOTO 140 DW(K) = I IF ( I .LE. N ) GOTO 130 IF ( K .GT. 1 ) GOTO 40 RETURN END C C SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DGESV computes the solution to a real system of linear equations * A * X = B, * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. * * The LU decomposition with partial pivoting and row interchanges is * used to factor A as * A = P * L * U, * where P is a permutation matrix, L is unit lower triangular, and U is * upper triangular. The factored form of A is then used to solve the * system of equations A * X = B. * * Arguments * ========= * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N coefficient matrix A. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension (N) * The pivot indices that define the permutation matrix P; * row i of the matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS matrix of right hand side matrix B. * On exit, if INFO = 0, the N-by-NRHS solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, so the solution could not be computed. * * ===================================================================== * * .. External Subroutines .. EXTERNAL DGETRF, DGETRS, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( NRHS.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGESV ', -INFO ) RETURN END IF * * Compute the LU factorization of A. * CALL DGETRF( N, N, A, LDA, IPIV, INFO ) IF( INFO.EQ.0 ) THEN * * Solve the system A*X = B, overwriting B with X. * CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, B, LDB, $ INFO ) END IF RETURN * * End of DGESV * END C C SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (preliminary version) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END C C SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. CHARACTER TRANS INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DGETRS solves a system of linear equations * A * X = B or A' * X = B * with a general N-by-N matrix A using the LU factorization computed * by DGETRF. * * Arguments * ========= * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A'* X = B (Transpose) * = 'C': A'* X = B (Conjugate transpose = Transpose) * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the right hand side matrix B. * On exit, the solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOTRAN * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLASWP, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NOTRAN = LSAME( TRANS, 'N' ) IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. $ LSAME( TRANS, 'C' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( NRHS.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRS', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * IF( NOTRAN ) THEN * * Solve A * X = B. * * Apply row interchanges to the right hand sides. * CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, 1 ) * * Solve L*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, $ ONE, A, LDA, B, LDB ) * * Solve U*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, $ NRHS, ONE, A, LDA, B, LDB ) ELSE * * Solve A' * X = B. * * Solve U'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, $ ONE, A, LDA, B, LDB ) * * Solve L'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, $ A, LDA, B, LDB ) * * Apply row interchanges to the solution vectors. * CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, -1 ) END IF * RETURN * * End of DGETRS * END C C SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DGETRF computes an LU factorization of a general M-by-N matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, IINFO, J, JB, NB * .. * .. External Subroutines .. EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRF', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN * * Use unblocked code. * CALL DGETF2( M, N, A, LDA, IPIV, INFO ) ELSE * * Use blocked code. * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Factor diagonal and subdiagonal blocks and test for exact * singularity. * CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) * * Adjust INFO and the pivot indices. * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO + J - 1 DO 10 I = J, MIN( M, J+JB-1 ) IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE * * Apply interchanges to columns 1:J-1. * CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) * IF( J+JB.LE.N ) THEN * * Apply interchanges to columns J+JB:N. * CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, $ IPIV, 1 ) * * Compute block row of U. * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) IF( J+JB.LE.M ) THEN * * Update trailing submatrix. * CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), $ LDA ) END IF END IF 20 CONTINUE END IF RETURN * * End of DGETRF * END C C SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END C C INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, $ N4 ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 * .. * * Purpose * ======= * * ILAENV is called from the LAPACK routines to choose problem-dependent * parameters for the local environment. See ISPEC for a description of * the parameters. * * This version provides a set of parameters which should give good, * but not optimal, performance on many of the currently available * computers. Users are encouraged to modify this subroutine to set * the tuning parameters for their particular machine using the option * and problem size information in the arguments. * * This routine will not function correctly if it is converted to all * lower case. Converting it to all upper case is allowed. * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies the parameter to be returned as the value of * ILAENV. * = 1: the optimal blocksize; if this value is 1, an unblocked * algorithm will give the best performance. * = 2: the minimum block size for which the block routine * should be used; if the usable block size is less than * this value, an unblocked routine should be used. * = 3: the crossover point (in a block routine, for N less * than this value, an unblocked routine should be used) * = 4: the number of shifts, used in the nonsymmetric * eigenvalue routines * = 5: the minimum column dimension for blocking to be used; * rectangular blocks must have dimension at least k by m, * where k is given by ILAENV(2,...) and m by ILAENV(5,...) * = 6: the crossover point for the SVD (when reducing an m by n * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds * this value, a QR factorization is used first to reduce * the matrix to a triangular form.) * = 7: the number of processors * = 8: the crossover point for the multishift QR and QZ methods * for nonsymmetric eigenvalue problems. * = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * =10: ieee NaN arithmetic can be trusted not to trap * =11: infinity arithmetic can be trusted not to trap * * NAME (input) CHARACTER*(*) * The name of the calling subroutine, in either upper case or * lower case. * * OPTS (input) CHARACTER*(*) * The character options to the subroutine NAME, concatenated * into a single character string. For example, UPLO = 'U', * TRANS = 'T', and DIAG = 'N' for a triangular routine would * be specified as OPTS = 'UTN'. * * N1 (input) INTEGER * N2 (input) INTEGER * N3 (input) INTEGER * N4 (input) INTEGER * Problem dimensions for the subroutine NAME; these may not all * be required. * * (ILAENV) (output) INTEGER * >= 0: the value of the parameter specified by ISPEC * < 0: if ILAENV = -k, the k-th argument had an illegal value. * * Further Details * =============== * * The following conventions have been used when calling ILAENV from the * LAPACK routines: * 1) OPTS is a concatenation of all of the character options to * subroutine NAME, in the same order that they appear in the * argument list for NAME, even if they are not used in determining * the value of the parameter specified by ISPEC. * 2) The problem dimensions N1, N2, N3, N4 are specified in the order * that they appear in the argument list for NAME. N1 is used * first, N2 second, and so on, and unused problem dimensions are * passed a value of -1. * 3) The parameter value returned by ILAENV is checked for validity in * the calling subroutine. For example, ILAENV is used to retrieve * the optimal blocksize for STRTRI as follows: * * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) * IF( NB.LE.1 ) NB = MAX( 1, N ) * * ===================================================================== * * .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER I, IC, IZ, NB, NBMIN, NX * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. External Functions .. INTEGER IEEECK EXTERNAL IEEECK * .. * .. Executable Statements .. * GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000, $ 1100 ) ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 100 CONTINUE * * Convert NAME to upper case if the first character is lower case. * ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN * * ASCII character set * IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF * ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN * * EBCDIC character set * IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) $ SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF * ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN * * Prime machines: ASCII+128 * IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF * C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) ) $ RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) * GO TO ( 110, 200, 300 ) ISPEC * 110 CONTINUE * * ISPEC = 1: block size * * In these examples, separate code is provided for setting NB for * real and complex. We assume that NB will take the same value in * single or double precision. * NB = 1 * IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN * 200 CONTINUE * * ISPEC = 2: minimum block size * NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 8 ELSE NBMIN = 8 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN * 300 CONTINUE * * ISPEC = 3: crossover point * NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN * 400 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 500 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 600 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 700 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 800 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * 900 CONTINUE * * ISPEC = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * ILAENV = 25 RETURN * 1000 CONTINUE * * ISPEC = 10: ieee NaN arithmetic can be trusted not to trap * ILAENV = 0 RETURN * 1100 CONTINUE * * ISPEC = 11: infinity arithmetic can be trusted not to trap * ILAENV = 0 RETURN * * End of ILAENV * END C C SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. INTEGER INCX, K1, K2, LDA, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DLASWP performs a series of row interchanges on the matrix A. * One row interchange is initiated for each of rows K1 through K2 of A. * * Arguments * ========= * * N (input) INTEGER * The number of columns of the matrix A. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the matrix of column dimension N to which the row * interchanges will be applied. * On exit, the permuted matrix. * * LDA (input) INTEGER * The leading dimension of the array A. * * K1 (input) INTEGER * The first element of IPIV for which a row interchange will * be done. * * K2 (input) INTEGER * The last element of IPIV for which a row interchange will * be done. * * IPIV (input) INTEGER array, dimension (M*abs(INCX)) * The vector of pivot indices. Only the elements in positions * K1 through K2 of IPIV are accessed. * IPIV(K) = L implies rows K and L are to be interchanged. * * INCX (input) INTEGER * The increment between successive values of IPIV. If IPIV * is negative, the pivots are applied in reverse order. * * Further Details * =============== * * Modified by * R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA * * ===================================================================== * * .. Local Scalars .. INTEGER I, I1, I2, INC, IP, IX, IX0, J, K, N32 DOUBLE PRECISION TEMP * .. * .. Executable Statements .. * * Interchange row I with row IPIV(I) for each of rows K1 through K2. * IF( INCX.GT.0 ) THEN IX0 = K1 I1 = K1 I2 = K2 INC = 1 ELSE IF( INCX.LT.0 ) THEN IX0 = 1 + ( 1-K2 )*INCX I1 = K2 I2 = K1 INC = -1 ELSE RETURN END IF * N32 = ( N / 32 )*32 IF( N32.NE.0 ) THEN DO 30 J = 1, N32, 32 IX = IX0 DO 20 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 10 K = J, J + 31 TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 10 CONTINUE END IF IX = IX + INCX 20 CONTINUE 30 CONTINUE END IF IF( N32.NE.N ) THEN N32 = N32 + 1 IX = IX0 DO 50 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 40 K = N32, N TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 40 CONTINUE END IF IX = IX + INCX 50 CONTINUE END IF * RETURN * * End of DLASWP * END C C SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DGETF2 computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the m by n matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER J, JP * .. * .. External Functions .. INTEGER IDAMAX EXTERNAL IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETF2', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * DO 10 J = 1, MIN( M, N ) * * Find pivot and test for singularity. * JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply the interchange to columns 1:N. * IF( JP.NE.J ) $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * INFO = J END IF * IF( J.LT.MIN( M, N ) ) THEN * * Update trailing submatrix. * CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, $ A( J+1, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of DGETF2 * END C C SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END C C LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END C C subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end C C subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end C C SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END C C integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end C C