* * $Id: sieval.F,v 1.1.1.1 1995/12/12 14:36:16 mclareni Exp $ * * $Log: sieval.F,v $ * Revision 1.1.1.1 1995/12/12 14:36:16 mclareni * Imported sources * * #include "sigma/pilot.h" *CMZ : 1.10/03 10/02/94 17.17.04 by Carlo E. Vandoni *-- Author : SUBROUTINE SIEVAL C .................................................. C PURPOSE C TO INTERPOLATE A FUNCTION FA=F(A) OVER A DIFFERENT SET OF C INDEPENDENT VARIABLE VALUES SO THAT RES=F(B) C C USAGE C SIGMA USER WRITES R=EVAL(FA,A,B) C C DESCRIPTION OF PARAMETERS C FA=F(A), THEIR NCO'S MUST BE IDENTICAL (MODES NEED NOT BE) C ROW LENGTH IS ARBITRARY C EXCEPT FOR ROW-LENGTH, NCO'S OF A AND B MUST BE COMPATIBLE C IN THE SENSE OF GENERALIZED ARITHMETIC C C RESULT HAS THE MODE OF FA THE ROW-LENGTH OF B C REST OF NCO DETERMINED BY GENERALIZED ARITHMETIC C C TRACKS USED C TRACK 38 FOR TRACING THE OUTER LOOP (ON FULL ROWS) C C COMM. BLOCKS USED C COM1 C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NGET C IGETT2 C SISTR2 C C METHOD C ROWS ARE COMBINED IN THE SENSE OF GENERALIZED ARITHMETIC C THE OUTER LOOP IS COPIED FROM OPER2 DO 50 IPC=... C C EACH POINT OF A ROW OF B IS COMPARED WITH C C FIRST THE LEFT END-POINT OF A ROW OF A WHERE C A MATCH IS OBTAINED WITHIN RELATIVE DIST. + OR - 1E-13 C C SECONDLY EACH INTERVAL (INCLUDING THE RIGHT END-POINT) C OF A ROW OF A C C NO MATCH PRODUCES A ZERO IN THE CORRESP. ELEM. OF RESULT C C A MATCH INTERPOLATES ON FA IN THE MATCHED INTERVAL C FOR COMPLEX A THE INTERVALS ARE RECTANGLES C IN THE COMPLEX PLANE C C REMARKS C ALL FOUR ARRAYS FA A B RES MUST BE IN CORE FOR GENERALIZD C ARITH LOOP TO WORK. THE RULE OF THREE IS NOT SUFFICIENT C ****************************************************************** C C C AUTHOR. JURIS REINFELDS DATE 10/02/75 C C C... PAW VERSION ... MAY 1988 C C .................................................. #include "sigma/sicsig.inc" #include "sigma/sigc.inc" #include "sigma/pawc.inc" C C C DIMENSION FNCO(10,3),IDOPE(11,3) DIMENSION B(2),ALEFT(2),ARIGHT(2) C DIMENSION IAD(3),JMIP(3),DIM(10) EQUIVALENCE(IAD(1),IZAD(1),IADA),(IADB,IZADB),(JMIP(1),IMODA,IZMOD ,(1)),(IMODB,IZMODB),(IMODC,IZMODC,MODE) C C CALL SITRAC(' SIEVAL ') C IN C=EVAL(FA,A,B) BOTH FA AND A ARE K=1, B IS K=2, C IS K=3 C C ICONCT DECIDES WHICH IS THE FIRST DIMENSION CONSIDERED C 1 FOR ALL DIMENSIONS C 2 FOR WHOLE ROWS ONLY C 3 FOR WHOLE MATRICES ONLY AND SO ON C C C IF ONE OR MORE INDICES MISSING DO 33 J=1,3 CALL SISTAK(J-1,MP,MN) IF(MN.EQ.MISIDX) GOTO 918 33 CONTINUE C C C IF THIRD ARG (B) IS NOT A NUMERICAL ITEM FNCO(1,2)=1.0 CALL SINGET(ISI,0,FNCO(1,2)) IF(IERRNO.NE.0)RETURN IF(ISI.GE.3) GOTO 9108 IDIMB=NDIM IZMODB=MODE IADB=IADDR-IZMODB C C C IF SECOND ARG (A) IS NOT A NUMERICAL ITEM FNCO(1,1)=1.0 CALL SINGET(ISI,1,FNCO(1,1)) IF(IERRNO.NE.0)RETURN IF(ISI.GE.3) GOTO 9108 IDIMA=NDIM IZMODA=MODE IADA=IADDR-IZMODA C MUST SAVE SCALAR VALUE IN CASE NEXT USE OF NGET IS ALSO A SCALAR XASCLR=DYNA(IADDR) XAIMAG=0. IF(IZMODA.EQ.2)XAIMAG=DYNA(IADDR+1) C C CHECK IF INCOMPATIBLE NCO S IN BINARY OPERATION C BUT OVERLOOK THE ROW-DIMENSION WHICH NEED NOT BE COMPATIBLE NDIM=MAX0(IDIMA,IDIMB) DO 2 I=2,NDIM IF((FNCO(I,1).EQ.FNCO(I,2)).OR.(FNCO(I,1).EQ.1.).OR.(FNCO(I,2).EQ. 11.)) GOTO 2 GOTO 968 2 CONTINUE C C IF FIRST ARG (FA) IS NOT A NUMERICAL ITEM DIM(1)=0.0 CALL SINGET(ISI,2,DIM) IF(IERRNO.NE.0)RETURN IF(ISI.GE.3) GOTO 9108 IADFA=IADDR-MODE C C CHECK IF FA AND A HAVE IDENTICAL NCO VECTORS IF(NDIM.NE.IDIMA) GOTO 968 DO 4 I=1,IDIMA IF(DYNA(IADFA+MODE-I).EQ.DYNA(IADA+IZMODA-I)) GOTO 4 GOTO 968 4 CONTINUE C C C MAKE NCO OF RESULT IN FNCO( ,3) C ROW-LENGTH OF RESULT IS ROWLENGTH OF B FNCO(1,3)=FNCO(1,2) ICONCT=2 C DO 20 I=ICONCT,NDIM FNCO(I,3)=AMAX1(FNCO(I,2),FNCO(I,1)) 20 CONTINUE C DO 24 I=NDIM,9 FNCO(I+1,3)=1. 24 CONTINUE C C MAKE ALL THREE DOPE VECTORS C MODE OF ANSWER FOLLOWS MODE OF FA C MODE IS EQUIVALENCED TO JMIP(3) DO 30 K=1,3 IDOPE(1,K)=JMIP(K) DO 130 J=1,10 IDOPE(J+1,K)=IDOPE(J,K)*IFIX(FNCO(J,K)) 130 CONTINUE 30 CONTINUE C C C GET AREA FOR THE RESULT LENGTH=IDOPE(NDIM+1,3) ISTRI=0 CALL SIGTT2(IADRES,LENGTH+NDIM,NDIM,FNCO(1,3)) IADDR=IADRES-MODE C C C NOW CALCULATE THE ROWLENGTH OF EACH ARRAY IROWA=IFIX(FNCO(1,1))*IZMODA-IZMODA IROWB=IFIX(FNCO(1,2))*IZMODB IROWC=IFIX(FNCO(1,3))*IZMODC B(2)=0. ALEFT(2)=0. ARIGHT(2)=0. C C NOW FILL EACH ROW OF ANSWER WITH ZEROS DO 40 I=MODE,LENGTH DYNA(IADDR+I)=0. 40 CONTINUE C C DO 50 IPC=IROWC,LENGTH,IROWC C MAIN LOOP FOR GENERALIZED ARRAY ARITHMETIC C COMBINES WHOLE ROWS (COMPATIBLE IN ALL DIMS EXCEPT ROW) I=ICONCT JRES=0 C DO 150 J=IZMODB,IROWB,IZMODB C TAKE EACH ELEMENT OF CURRENT ROW OF B JRES=JRES+MODE B(1)=DYNA(IADB+J) IF(IZMODB.EQ.2) B(2)=DYNA(IADB+J+1) C DO 250 K=IZMODA,IROWA,IZMODA C SEARCH EACH INTERVAL OF A INCLUDING RIGHT END POINT C BUT EXCLUDING LEFT END POINT ALEFT(1)=DYNA(IADA+K) ARIGHT(1)=DYNA(IADA+K+IZMODA) IF(IZMODA.NE.2) GOTO 251 ALEFT(2)=DYNA(IADA+K+1) ARIGHT(2)=DYNA(IADA+K+3) 251 CONTINUE C IF(K.NE.IZMODA) GOTO 2518 C NOW HANDLE THE LEFT-END-POINT AS A SPECIAL CASE C END-POINT EQUALITY COUNTS WITHIN RELATIVE DISTANCE OF 1+OR-1E-13 C FOR SCALAR A AND FA, DYNA(IADA+IZMODA) CONTAINS FA NOT A IF(IDIMA.NE.1.OR.IROWA.GT.0) GOTO 2511 ALEFT(1)=XASCLR ALEFT(2)=XAIMAG 2511 CONTINUE C DO 2512 KK=1,2 IF(B(KK).EQ.ALEFT(KK))GOTO 2512 IF(ABS(B(KK)-ALEFT(KK))/AMAX1(ABS(B(KK)),ABS(ALEFT(KK))) ,.GT.1.E-13) GOTO 2516 2512 CONTINUE C C SET LEFT END POINT INTO ANSWER AND SEARCH NO FURTHER KRES=IADFA+K*MODE/IZMODA DYNA(IADDR+JRES)=DYNA(KRES) IF(MODE.EQ.2) DYNA(IADDR+JRES+1)=DYNA(KRES+1) GOTO 150 C 2516 CONTINUE C SKIP INTERVAL PROCESSING FOR ROW-LENGTH ONE IF(IROWA.LE.0)GOTO 150 C 2518 CONTINUE C DO 252 KK=1,2 C COMPARE REAL AND IMAGINARY PARTS OF B-POINT WITH INTERVAL OF A C IF BOTH END POINTS ARE EQUAL TO B-POINT, WE HAVE A MATCH IF(B(KK).EQ.ALEFT(KK).AND.B(KK).EQ.ARIGHT(KK)) GOTO 252 C IF POINT IS OUTSIDE THE INTERVAL, WE HAVE NO MATCH, TRY NEXT INTER IF((B(KK).LE.ALEFT(KK).OR.B(KK).GT.ARIGHT(KK)).AND. ,(B(KK).GE.ALEFT(KK).OR.B(KK).LT.ARIGHT(KK))) GOTO 250 252 CONTINUE C C GET HERE ONLY IF POINT IS WITHIN INTERVAL SO TRANSFER FA TO RESULT KRES=K*MODE/IZMODA DO 254 KK=1,2 DYNA(IADDR+JRES+KK-1)=DYNA(IADFA+KRES+KK-1) C IF POINT NOT IN INTERVAL WHOSE END POINTS ARE EQUAL IF(B(KK).NE.ALEFT(KK)) DYNA(IADDR+JRES+KK-1)= ,(DYNA(IADFA+KRES+KK-1)-DYNA(IADFA+KRES+MODE+KK-1))* ,(B(KK)-ARIGHT(KK))/(ALEFT(KK)-ARIGHT(KK))+DYNA(IADFA+KRES+MODE ,+KK-1) C IF POINT IS IN ONE INTERVAL, LOOK NO FURTHER, GET NEXT POINT IF(MODE.NE.2) GOTO 150 254 CONTINUE GOTO 150 C 250 CONTINUE C 150 CONTINUE IF(SITRAK(38).NE.1) GOTO 154 C FIRST ITEM IS THE LAST ELEMENT OF PREVIOUS ROW IENDFA=IADFA+(IROWA+IZMODA)*MODE/IZMODA IENDA=IADA+IROWA+IZMODA IENDB=IADB+IROWB IENDC=IADDR+IROWC WRITE(NPRINT,1001) (DYNA(I),I=IADFA,IENDFA) WRITE(NPRINT,1002) (DYNA(I),I=IADA,IENDA) WRITE(NPRINT,1003) (DYNA(I),I=IADB,IENDB) WRITE(NPRINT,1004) (DYNA(I),I=IADDR,IENDC) 1001 FORMAT(' FA=',10F6.1) 1002 FORMAT(' A=',10F6.1) 1003 FORMAT(' B=',10F6.1) 1004 FORMAT(' C=',10F6.1) C 154 CONTINUE C C DO 156 J=1,3 C STEP THE ADDRESSES FORWARD BY ONE ROW IAD(J)=IAD(J)+IFIX(FNCO(1,J))*JMIP(J) 156 CONTINUE IADFA=IADFA+IFIX(FNCO(1,1))*MODE C 55 CONTINUE C CHECK WHICH SUBSTRUCTURE IS COMPLETED IF(MOD(IPC,IDOPE(I+1,3)).NE.0) GOTO 57 I=I+1 IF(IPC.LT.LENGTH) GOTO 55 C C 57 CONTINUE DO 157 K=1,2 C REPEAT A SUBSTRUCTURE IF NECESSARY IF(FNCO(I,K).EQ.1.) IAD(K)=IAD(K)-IDOPE(I,K) 157 CONTINUE IF(FNCO(I,1).EQ.1.) IADFA=IADFA-IDOPE(I,1)*MODE/IZMODA C 50 CONTINUE C C C STORE REQUIRES THAT IADDR POINTS TO REAL PART OF ONE-COMP ARRAY IADDR=IADRES CALL SISTR2(3) RETURN C 918 CONTINUE C MISSING INDEX IS MEANINGLESS IN EVAL CNAME='EVAL ' CALL SINERR(18) RETURN C 968 CONTINUE C INCOMPATIBLE ARGUMENTS CNAME='EVAL ' CALL SINERR(68) RETURN C 9108 CONTINUE C ONE ARG IS A PROGRAM NAME CNAME='EVAL ' CALL SINERR(58) C C 999 END