C
SUBROUTINE SOFT * (NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,WEIGHT,NM,INDM,CS,MB,B) C INTEGER NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,NM,INDM(*),MB REAL WEIGHT,CS(*),B(MB) C C This subroutine accumulates the prior subjective information C covariance matrix describing the smoothness of the functions C interpolated by means of subroutines of the file 'val.for'. C C Input: C NCLASS..Number of classes required. C NA1,NA2,NA3... Orders of X1,X2,X3-partial derivatives of the first C basis function in the product being integrated. C NB1,NB2,NB3... Orders of X1,X2,X3-partial derivatives of the C second basis function in the product being integrated. C WEIGHT..Weighting factor corresponding to the derivatives. C CS... Symmetric matrix to be regularized: C Contribution corresponding to the derivatives will be C added to the matrix CS if WEIGHT.NE.0. C Not used if WEIGHT=0. C MB... Dimension of working array B. C B... Working array. Auxiliary storage locations for partial C covariance matrices. For the comments on its minimum C dimension refer to the destription of error C 372. Not used if WEIGHT=0. C C Output: C NM... Number of the unknown model parameters. C INDM... Indices of the unknown model parameters. Integer array of C dimension NM. C CS... Resulting symmetric subjective prior information C covariance matrix (input matrix increased by WEIGHT times C the product of the derivatives averaged over the B-spline C grid). Real array of dimension NM*(NM+1)/2. C Unchanged input if WEIGHT=0. C C Common block: INCLUDE 'val.inc' C val.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SPSP C C Date: 1997, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER ICLASS,JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,NVAR INTEGER NX(3),NX1,NX2,NX3 EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3) INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7 EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3) EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6) EQUIVALENCE (JADR(7),JADR7) INTEGER I1,I2,I3,II,I,J,N REAL WDENS,SIGMA C C B... Auxiliary storage locations for partial covariance C matrices. C ICLASS..Index of the class. C JGROUP..Address of a current group. C LFUNCT,MFUNCT,JFUNCT... Addresses of the first, last and arbitrary C functions of the group. C LADR,MADR,IADR... Addresses of the first, last and arbitrary C parameters of the current function. C NVAR... Number of the independent variables A1, A2, A3 of the C interpolated function W. C NX=(NX1,NX2,NX3)... Numbers of grid lines. C JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... Addresses of C parameters describing the interpolated function (grid C coordinates, B-spline coefficients, B-spline basis C functions). C I1,I2,I3,I,J,N... Local auxiliary variables. C WDENS...Weighting factor divided by volume (area,length). C SIGMA...Tension factor. C C....................................................................... C NM=0 IF(NCLASS.LT.1.OR.IPAR(0).LT.NCLASS) THEN C 371 PAUSE 'Error 371 in SOFT: Incorrect number of classes' STOP C The number of classes required is zero, negative or greater than C the number of classes defined. END IF C Loop for classes: DO 92 ICLASS=1,NCLASS C Loop for groups of the current class: DO 91 JGROUP=IPAR(ICLASS-1)+1,IPAR(ICLASS) LFUNCT=IPAR(JGROUP-1)+1 MFUNCT=IPAR(JGROUP) MADR =IPAR(LFUNCT-1) C C Loop for functions of the current group: DO 90 JFUNCT=LFUNCT,MFUNCT C Starting and end addresses of the parameters describing the C function LADR=MADR+1 MADR=IPAR(JFUNCT) IF(LADR.LE.MADR) THEN WDENS=WEIGHT C Tension factor SIGMA=RPAR(LADR+5) C C The number, types and values of the independent variables C of the interpolated function: C Initial address IADR=LADR+6 C Initial number of the independent variables NVAR=0 NX1=1 NX2=1 NX3=1 JADR1=0 JADR2=0 JADR3=0 JADR4=0 C Loop for the possible independent variables: DO 20 I=LADR+2,LADR+4 C Type of the possible independent variable: J=IPAR(I) IF(J.GT.0) THEN N=IABS(IPAR(IADR)) IF(N.GE.2) THEN NVAR=NVAR+1 NX(NVAR)=N ELSE IF(N.EQ.1) THEN JADR(NVAR+1)=JADR(NVAR+1)+1 END IF IADR=IADR+1 END IF 20 CONTINUE JADR5=0 JADR6=0 JADR7=0 C C Interpolated function W: JADR1=IADR+JADR1 IF(NVAR.LE.0) THEN C No independent variable. JADR4=JADR1 ELSE JADR2=JADR1+NX1+JADR2 WDENS=WDENS/(RPAR(JADR1+NX1-1)-RPAR(JADR1)) IF(NVAR.EQ.1) THEN C One independent variable: JADR4=JADR2 JADR5=JADR2+NX1 ELSE JADR3=JADR2+NX2+JADR3 WDENS=WDENS/(RPAR(JADR2+NX2-1)-RPAR(JADR2)) IF(NVAR.EQ.2) THEN C Two independent variables: JADR4=JADR3 JADR5=JADR3+NX1*NX2 JADR6=JADR4+5*NX1 ELSE C Three independent variables: JADR4=JADR3+NX3+JADR4 JADR5=JADR4+NX1*NX2*NX3 JADR6=JADR5+5*NX1 JADR7=JADR6+5*NX2 WDENS=WDENS/(RPAR(JADR3+NX3-1)-RPAR(JADR3)) END IF END IF END IF C WDENS is weighting factor divided by volume (area,length). C JADR4=JADR4-1 N=NX1*NX2*NX3 IF(WEIGHT.EQ.0.) THEN DO 72 I2=1,N INDM(NM+I2)=JADR4+I2 72 CONTINUE NM=NM+N ELSE IF(NA1.EQ.NB1.AND.NA2.EQ.NB2.AND.NA3.EQ.NB3) THEN I=N*(N+1)/2 J=0 ELSE I=N*N J=1 END IF IF(NA1.EQ.NB1) THEN I1=NX1*(NX1+1)/2 ELSE I1=NX1*NX1 END IF IF(NA2.EQ.NB2) THEN I2=NX2*(NX2+1)/2 ELSE I2=NX2*NX2 END IF IF(NA3.EQ.NB3) THEN I3=NX3*(NX3+1)/2 ELSE I3=NX3*NX3 END IF I1=I1+I I2=I2+I1 I3=I3+I2 IF(I3.GT.MB) THEN C 372 PAUSE 'Error 372 in SOFT: Insufficient working memory' STOP C The dimension MB of the buffer B is not sufficient, C see the above command lines, where NX1,NX2,NX3 are C the dimensions of the grid for spline interpolation. C MB and B(MB) are the dummy arguments of this C subroutine. If the actual argument is allocated in C array RAM(MRAM) declared in include file C ram.inc, C you may wish to increase MRAM. END IF DO 80 II=1,I B(II)=0. 80 CONTINUE CALL SPSP(NA1,NA2,NA3,NB1,NB2,NB3,NX1,NX2,NX3, * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3), * RPAR(JADR5),RPAR(JADR6),RPAR(JADR7), * SIGMA,WDENS,B(1),B(I+1),B(I1+1),B(I2+1)) I=NM*(NM+1)/2 DO 82 I2=1,N INDM(NM+I2)=JADR4+I2 I=I+NM DO 81 I1=1,I2 I=I+1 IF(J.EQ.0) THEN CS(I)=CS(I)+B(I2*(I2-1)/2+I1) ELSE CS(I)=CS(I)+B(N*(I2-1)+I1)+B(N*(I1-1)+I2) END IF 81 CONTINUE 82 CONTINUE NM=NM+N C Contribution corresponding to the derivatives C is added to the matrix CS. END IF C END IF 90 CONTINUE 91 CONTINUE 92 CONTINUE C End of loops for functions. C RETURN END C C======================================================================= C