C
SUBROUTINE SOFT * (ICLASS,NA1,NA2,NA3,NB1,NB2,NB3,NW,WEIGHT,NM,INDM,CS,MB,B) C INTEGER ICLASS,NA1,NA2,NA3,NB1,NB2,NB3,NW,NM,INDM(*),MB REAL WEIGHT(NW,ICLASS),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 ICLASS..Index of the class 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 If triplets (NA1,NA2,NA3) and (NB1,NB2,NB3) are C different, both Sobolev scalar products corresponding to C derivatives (NA1,NA2,NA3),(NB1,NB2,NB3) and C (NA1,NA2,NA3),(NB1,NB2,NB3) are added to matrix CS. C NW... Maximum index that identifies the physical meaning of the C function. C NM... Input number of the unknown model parameters. C Zero for the first required class, output from the C previous invocation for subsequent classes of functions. C WEIGHT..Weighting factors corresponding to the derivatives. C WEIGHT(MPAR,ICLASS) is the weighting factor corresponding C to the function indexed MPAR of the ICLASSth class. C CS... Symmetric matrix to be regularized: C Contribution corresponding to the derivatives will be C added to the corresponding part of matrix CS if some of C WEIGHT(MPAR,ICLASS).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... Output 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). If triplets (NA1,NA2,NA3) and (NB1,NB2,NB3) are C different, both Sobolev scalar products corresponding to C derivatives (NA1,NA2,NA3),(NB1,NB2,NB3) and C (NA1,NA2,NA3),(NB1,NB2,NB3) are added. C CS is stored as a 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: 2000, July 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,NVAR INTEGER NA(3),NB(3),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,MPAR 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,MPAR... Local auxiliary variables. C WDENS...Weighting factor divided by volume (area,length). C SIGMA...Tension factor. C C....................................................................... C C NM=0 C IF(NCLASS.LT.1.OR.IPAR(0).LT.NCLASS) THEN IF(ICLASS.LT.1.OR.IPAR(0).LT.ICLASS) THEN C 371 CALL ERROR('371 in SOFT: Incorrect number of classes') C The number of classes required is zero, negative or greater than C the number of classes defined. END IF C Loop for classes: C 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 MPAR=IPAR(LADR) IF(MPAR.LT.1.OR.MPAR.GT.NW) THEN C 373 CALL ERROR('373 in SOFT: Incorrect index of function') C The index that identifies the physical meaning of the C function is not in range 1 to MPAR. Check the C invocations of subroutine VAL1 and of this subroutine. END IF WDENS=WEIGHT(MPAR,ICLASS) 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 NA(1)=0 NA(2)=0 NA(3)=0 NB(1)=0 NB(2)=0 NB(3)=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 IF(J.EQ.1) THEN NA(NVAR)=NA1 NB(NVAR)=NB1 ELSE IF(J.EQ.2) THEN NA(NVAR)=NA2 NB(NVAR)=NB2 ELSE IF(J.EQ.3) THEN NA(NVAR)=NA3 NB(NVAR)=NB3 END IF 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=JADR4+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=JADR4+NX1*NX2 JADR6=JADR5+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 I=NA1+NA2+NA3-NA(1)-NA(2)-NA(3) * +NB1+NB2+NB3-NB(1)-NB(2)-NB(3) IF(WDENS.EQ.0..OR.I.NE.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(NA(1).EQ.NB(1)) THEN I1=NX1*(NX1+1)/2 ELSE I1=NX1*NX1 END IF IF(NA(2).EQ.NB(2)) THEN I2=NX2*(NX2+1)/2 ELSE I2=NX2*NX2 END IF IF(NA(3).EQ.NB(3)) 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 CALL ERROR('372 in SOFT: Insufficient working memory') 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 of a single function describing a surface or a C material parameter. Dimension MB must be sufficiently C large for each interpolated function. 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 C ram.inc, you may wish to increase MRAM. END IF DO 80 II=1,I B(II)=0. 80 CONTINUE CALL SPSP(NA(1),NA(2),NA(3),NB(1),NB(2),NB(3), * 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 C 92 CONTINUE C End of loops for functions. C RETURN END C C======================================================================= C