1 : C LEVMAR SOURCE FANDEUR 10/07/07 21:15:17 6708 2 : SUBROUTINE LEVMAR 3 : * 4 : * moteur pour algorithme levenberg-marquardt 5 : * 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8(A-H,O-Z) 8 : 9 : -INC CCOPTIO 10 : -INC SMLENTI 11 : -INC SMLREEL 12 : SEGMENT MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL,MLREE8.MLREEL, 13 : + MLREE9.MLREEL,MLRE11.MLREEL,MLRE12.MLREEL,MLRE13.MLREEL 14 : 15 : PARAMETER (NTYP=11) 16 : CHARACTER*4 LISTYP(NTYP) 17 : DATA LISTYP / 'ABSC', 'ORDO', 'PARA', 'TORD', 'DYDA', 18 : & 'SIGM', 'PROC', 'ALPH', 'BETA', 'CHI2', 19 : & 'PARO' / 20 : 21 : data iprocc/0/ 22 : SAVE iprocc 23 : c 24 : write(ioimp,*) 'levm ', iprocc 25 : c 26 : jg = NTYP 27 : segini mlenti 28 : c 29 : if (iprocc.gt.0) then 29 : call LIROBJ('LISTREEL',ip1,1,iret) 30 : IF (IERR.NE.0) GOTO 900 31 : lect(4) = ip1 31 : call LIROBJ('LISTREEL',ip1,1,iret) 32 : IF (IERR.NE.0) GOTO 900 33 : lect(5) = ip1 34 : endif 35 : c 36 : idon = 0 37 : 10 continue 37 : CALL LIRMOT(LISTYP,NTYP,iplac,0) 38 : IF (IERR.NE.0) GOTO 900 39 : if (iplac.gt.0) then 40 : if (iplac.ne.7) then 40 : call LIROBJ('LISTREEL',ip1,1,iret) 41 : else if (iplac.eq.7) then 41 : call LIROBJ('PROCEDUR',ip1,1,iret) 42 : endif 43 : IF (IERR.NE.0) GOTO 900 44 : lect(iplac) = ip1 45 : idon = idon + 1 46 : goto 10 47 : endif 48 : 49 : if (iprocc.eq.0) then 50 : if (lect(7).le.0) then 51 : c* iprocc = 0 52 : moterr(1:8) ='PROCEDUR' 52 : call erreur(37) 53 : GOTO 900 54 : endif 55 : mlree1 = lect(1) 56 : mlree2 = lect(2) 57 : mlree6 = lect(6) 58 : mlree3 = lect(3) 59 : if (mlree1.le.0 .or. mlree2.le.0 .or. mlree6.le.0 .or. 60 : & mlree3.le.0) then 61 : c* iprocc = 0 61 : call erreur(201) 62 : GOTO 900 63 : endif 64 : segact,mlree1,mlree2,mlree6 65 : n1 = mlree1.prog(/1) 66 : n2 = mlree2.prog(/1) 67 : n6 = mlree6.prog(/1) 68 : if (n1.LE.0 .or. n1.NE.n2 .or. n1.NE.n6 .or. n2.NE.n6) then 69 : segdes,mlree1,mlree2,mlree6 69 : call ERREUR(217) 70 : c* iprocc = 0 71 : GOTO 900 72 : endif 73 : segact,mlree3 74 : ma = mlree3.prog(/1) 75 : segdes,mlree3 76 : if (ma.LE.0) then 77 : segdes,mlree1,mlree2,mlree6 77 : call ERREUR(21) 78 : GOTO 900 79 : endif 80 : segini,mlre11=mlree3 81 : lect(3) = mlre11 82 : jg = ma 83 : segini mlre12 84 : lect(11) = mlre12 85 : jg = ma*ma 86 : segini mlree8 87 : lect(8) = mlree8 88 : jg = ma 89 : segini mlree9 90 : lect(9) = mlree9 91 : jg = 2 92 : segini mlreel 93 : prog(1) = -1. 94 : prog(2) = -1. 95 : lect(10) = mlreel 96 : c 97 : else 98 : c* else if (iprocc.gt.0) then 99 : do idon = 1, NTYP 100 : if (lect(idon).LE.0) then 101 : write(ioimp,*) '!! IPROCC =',iprocc 102 : iprocc = 0 103 : if (idon.eq.7) then 104 : moterr(1:8) ='PROCEDUR' 104 : call erreur(37) 105 : else 105 : call erreur(201) 106 : endif 107 : GOTO 900 108 : endif 109 : enddo 109 : call levma1(mlenti) 110 : endif 111 : 112 : mlreel = lect(10) 113 : chi2 = prog(1) 114 : alamda = prog(2) 115 : * 116 : if(iprocc.gt.100.or.abs(chi2).lt.1.e-5.or.alamda.gt.1.e11) goto 80 117 : * 118 : iprocc = iprocc + 1 119 : ip1 = lect(11) 119 : call nomobj('LISTREEL','WLEVPARO',ip1) 119 : call ecrobj('LISTREEL',ip1) 119 : call ECRCHA('PARO') 120 : ip1 = lect(9) 120 : call nomobj('LISTREEL','WLEVBETA',ip1) 120 : call ecrobj('LISTREEL',ip1) 120 : call ECRCHA('BETA') 121 : ip1 = lect(8) 121 : call nomobj('LISTREEL','WLEVALPH',ip1) 121 : call ecrobj('LISTREEL',ip1) 121 : call ECRCHA('ALPH') 122 : ip1 = lect(7) 122 : call ecrobj('PROCEDUR',ip1) 122 : call ECRCHA('PROC') 123 : ip1 = lect(10) 123 : call nomobj('LISTREEL','WLEVCHI2',ip1) 123 : call ecrobj('LISTREEL',ip1) 123 : call ECRCHA('CHI2') 124 : ip1 = lect(6) 124 : call ecrobj('LISTREEL',ip1) 124 : call ECRCHA('SIGM') 125 : ip1 = lect(3) 125 : call nomobj('LISTREEL','WLEVPARA',ip1) 125 : call ecrobj('LISTREEL',ip1) 125 : call ECRCHA('PARA') 126 : ip1 = lect(2) 126 : call ecrobj('LISTREEL',ip1) 126 : call ECRCHA('ORDO') 127 : ip1 = lect(1) 127 : call ecrobj('LISTREEL',ip1) 127 : call ECRCHA('ABSC') 127 : call ECRCHA('LEVM') 128 : ip1 = lect(3) 128 : call ecrobj('LISTREEL',ip1) 129 : ip1 = lect(1) 129 : call ecrobj('LISTREEL',ip1) 130 : ip1 = lect(7) 130 : call ecrobj('PROCEDUR',ip1) 131 : GOTO 900 132 : 133 : 80 continue 134 : if (alamda.gt.1.e11) write(ioimp,*) 'je ne sais pas faire mieux' 135 : * oubli des objets intermediares 136 : ip1 = lect(8) 136 : CALL NOMOBJ('ANNULE ','WLEVALPH',ip1) 137 : ip1 = lect(9) 137 : CALL NOMOBJ('ANNULE ','WLEVBETA',ip1) 138 : ip1 = lect(11) 138 : CALL NOMOBJ('ANNULE ','WLEVPARO',ip1) 139 : 140 : mlree1 = lect(1) 141 : mlree2 = lect(2) 142 : mlree3 = lect(6) 143 : segdes mlree1,mlree2,mlree3 144 : mlree1 = lect(3) 145 : segini,mlree2=mlree1 146 : ip1 = mlree1 146 : CALL NOMOBJ('ANNULE ','WLEVPARA',ip1) 147 : ipout = mlree2 148 : mlreel = lect(10) 149 : chi2 = prog(1) 150 : ip1 = mlreel 150 : CALL NOMOBJ('ANNULE ','WLEVPARA',ip1) 150 : call ECRREE(chi2) 150 : call ECROBJ('LISTREEL', ipout) 151 : iprocc = 0 152 : 153 : 900 CONTINUE 154 : segsup mlenti 155 : 156 : RETURN 157 : END 158 : 159 :
© Cast3M 2003 - All rights reserved.
Disclaimer