OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
lag_mult_solv.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "lagmult.inc"
#include "com08_c.inc"
#include "units_c.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine lag_mult_solv (nh, nc, ncr, a, v, mas, iadll, lll, jll, xll, iadh, jcih, hh, z, p, r, q, ltsm, hl, diag_h, diag_l, work1, work2, work3, lambda, rbyl, npbyl, ar, vr, in, iadhf, jcihf, icftag, jcftag, ncf_s, ncf_e)
subroutine lag_mult_solvp (nh, nc, ncr, a, v, mas, iadll, lll, jll, xll, iadh, jcih, hh, z, p, r, q, ltsm, hl, diag_h, diag_l, work1, work2, work3, lambda, rbyl, npbyl, ar, vr, in, iadhf, jcihf, icftag, jcftag, ncf_s, ncf_e, indexlag)
subroutine lag_mult_sdp (nh, nc, ncr, a, v, mas, iadll, lll, jll, xll, iadh, jcih, hh, z, p, r, q, ltsm, hl, diag_h, diag_l, work1, work2, work3, lambda, rbyl, npbyl, ar, vr, in, iadhf, jcihf, icftag, jcftag, ncf_s, ncf_e, indexlag)

Function/Subroutine Documentation

◆ lag_mult_sdp()

subroutine lag_mult_sdp ( integer nh,
integer nc,
integer ncr,
a,
v,
mas,
integer, dimension(*) iadll,
integer, dimension(*) lll,
integer, dimension(*) jll,
xll,
integer, dimension(*) iadh,
integer, dimension(*) jcih,
hh,
z,
p,
r,
q,
ltsm,
hl,
diag_h,
diag_l,
work1,
integer, dimension(*) work2,
integer, dimension(*) work3,
lambda,
rbyl,
integer, dimension(nnpby,*) npbyl,
ar,
vr,
in,
integer, dimension(*) iadhf,
integer, dimension(*) jcihf,
integer, dimension(*) icftag,
integer, dimension(*) jcftag,
integer ncf_s,
integer ncf_e,
integer, dimension(*) indexlag )

Definition at line 793 of file lag_mult_solv.F.

802C-----------------------------------------------
803C I m p l i c i t T y p e s
804C-----------------------------------------------
805#include "implicit_f.inc"
806C-----------------------------------------------
807C C o m m o n B l o c k s
808C-----------------------------------------------
809#include "param_c.inc"
810#include "task_c.inc"
811#include "com08_c.inc"
812#ifdef MUMPS5
813#include "dmumps_struc.h"
814#include "lagmult.inc"
815#endif
816C-----------------------------------------------
817C D u m m y A r g u m e n t s
818C-----------------------------------------------
819 INTEGER NC,NCR,NCF_S,NCF_E,NH
820 INTEGER LLL(*),JLL(*),IADLL(*),IADH(*),JCIH(*),IADHF(*),
821 . JCIHF(*),NPBYL(NNPBY,*),WORK2(*),WORK3(*),
822 . ICFTAG(*),JCFTAG(*), INDEXLAG(*)
823 my_real
824 . mas(*),in(*), a(3,*),ar(3,*),v(3,*),vr(3,*),
825 . xll(*),p(*),r(*),q(*),lambda(*),rbyl(nrby,*),
826 . ltsm(6,*),z(*),hh(*),hl(*),diag_h(*),diag_l(*),
827 . work1(*)
828C-----------------------------------------------
829C L o c a l V a r i a b l e s
830C-----------------------------------------------
831 INTEGER I,J,IK,IC,JC,IH,IP,ITER,NITER,ITERMAX,NIN,LENH,NZ
832 my_real
833 . as,asmax,alpha,beta,r2,r2new,rnorm,pq,hij,scale,dd
834#ifdef MUMPS5
835 type(dmumps_struc)::MUMPS_PAR
836#endif
837C-----------------------------------------------
838C External Functions
839C-----------------------------------------------
840 INTEGER CHOLFACT
841 EXTERNAL cholfact
842C-----------------------------------------------
843C
844C | M | Lt| | a | M ao
845C |---+---| | = |
846C | L | 0 | | la | bo
847C
848C [M] a + [L]t la = [M] ao
849C [L] a = bo
850C
851C a = -[M]-1[L]t la + ao
852C [L][M]-1[L]t la = [L] ao - bo
853C
854C [H] = [L][M]-1[L]t
855C b = [L] ao - bo
856C [H] la = b
857C
858C a = ao - [M]-1[L]t la
859C-----------------------------------------------
860C
861C la : LAMBDA(NC)
862C ao : A(3,NUMNOD)
863C [M](i) : MAS(NUMNOD)
864C [L][M]-1[L]t la : Q(NC)
865C [L] ao - b : B(NC)
866C [M]-1[L]t : LTSM(3,NUMNOD) une colonne de [M]-1[L]t
867C : Z(NC)
868C
869C [L](ic,j,i) -> XLL(IK)
870C ic = 1..nc : IK = IADLL(IC)...IADLL(IC+1)-1
871C j = 1,2,3 : J = JLL(IK)
872C i = 1..numnod : I = LLL(IK)
873C
874C
875C NC : nombre de contact
876C
877C IC : numero du contact (1,NC)
878C IK : index de XLL,LLL,JLL
879C I : numero global du noeud (1,NUMNOD)
880C J : direction 1,2,3
881C
882C-----------------------------------------------
883C [H](ic,jc) -> HH(IH)
884C ic = 1..nc : IH = IADH(IC)...IADH(IC+1)-1
885C jc = 1..ic : JC = JCIH(IH)
886C
887C q = [H] p
888C
889C do ic=1,nc
890C q(ic) = 0
891C do jc=1,nc
892C q(ic) = q(ic) + h(ic,jc)*p(jc)
893C enddo
894C enddo
895C
896C DO IC=1,NC
897C Q(IC) = 0
898C DO IH=IADH(IC),IADH(IC+1)-1
899C Q(IC) = Q(IC) - HH(IH)*P(JCIH(IH))
900C ENDDO
901C ENDDO
902C-----------------------------------------------
903C evaluation de b:
904C
905C Vs = Somme(Ni Vi) + vout
906C Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai) + vout
907C Somme(dt Ni Ai) - dt As = Vs_ -Somme(Ni Vi_) - vout
908C [L] = dt {N1,N2,..,N15,-1}
909C bo = {N1,N2,..,N15,-1} a = -[L]/dt v_ - vout
910C b = [L] ao - bo
911C b = [L] ao + [L]/dt v_ + vout
912C b = [L]/dt vo+ + vout
913C-----------------------------------------------------------------------
914C
915C calcul de [H] la = b par gradient conjugue
916C
917C-----------------------------------------------------------------------
918C Gradient conjugue [H] la = b avec preconditionnement [N]
919C [H]-1 = Somme([I-H]^i),i=0,1...inf si |I-H|<1
920C [N]-1 = [2I-H]
921C-----------------------------------------------------------------------
922C r(0) = b - [H] la(0)
923C for i=1,niter
924C z(i-1) = [N]-1 r(i-1)
925C rho(i-1) = r(i-1)T z(i-1)
926C if (i == 1)
927C p(1) = z(0)
928C else
929C beta = rho(i-1) / rho(i-2)
930C p(i) = z(i-1) + beta p(i-1)
931C endif
932C q(i) = [H] p(i)
933C alpha = rho(i-1) / p(i)T q(i)
934C la(i) = la(i-1) + alpha p(i)
935C r(i) = r(i-1) - alpha q(i)
936C end
937C======================================================================|
938#ifdef MUMPS5
939 IF(ispmd==0) THEN
940 scale = three_over_4
941 asmax = zep6
942 r2new = zero
943 lag_nc= nc
944C---
945C--- Calcul de la matrice H
946 CALL lag_mult_hp(nc ,lenh ,nh ,mas ,in ,
947 2 diag_h ,hh ,iadll ,lll ,jll ,
948 3 xll ,ltsm ,iadhf ,jcihf ,iadh ,
949 4 jcih ,rbyl ,npbyl ,icftag ,jcftag ,
950 5 ncf_s ,ncf_e ,ncr ,indexlag)
951C------------------------------------------
952C Normalisation du systeme
953c lagopt=0 => sans norme
954c lagopt=1 => norme L1 (diagonal)
955c lagopt=2 => norme L2 (columns)
956C------------------------------------------
957 IF (lagopt==1) THEN
958 dd = zero
959 DO ic=1,nc
960 dd = dd + diag_h(ic)*diag_h(ic)
961 ENDDO
962 DO ic=1,nc
963 z(ic) = dd
964 ENDDO
965 ELSEIF (lagopt==2) THEN
966 DO ic=1,nc
967 z(ic) = diag_h(ic)*diag_h(ic)
968 ENDDO
969 DO ic=1,nc
970 DO ih=iadh(ic),iadh(ic+1)-1
971 jc = jcih(ih)
972 hij = hh(ih)*hh(ih)
973 z(ic) = z(ic) + hij
974 z(jc) = z(jc) + hij
975 ENDDO
976 ENDDO
977 ENDIF
978C---
979 IF (lagopt>0) THEN
980 DO ic=1,nc
981 z(ic) = scale/sqrt(sqrt(z(ic)))
982 ENDDO
983C--- normalisation [H] Somme([H](ic,i)^2)=1 i=1...nc
984 DO ic=1,nc
985 diag_h(ic) = diag_h(ic)*z(ic)*z(ic)
986 DO ih=iadh(ic),iadh(ic+1)-1
987 hh(ih) = hh(ih)*z(ic)*z(jcih(ih))
988 ENDDO
989 ENDDO
990C--- normalisation [LLL] et {R}
991 DO ic=1,nc
992 DO ik=iadll(ic),iadll(ic+1)-1
993 xll(ik) = xll(ik)*z(ic)
994 ENDDO
995 r(ic) = r(ic)*z(ic)
996 ENDDO
997 ENDIF
998C------------------------------------------
999C Factorisation de [H] : Cholesky incomplet + shift iteratif
1000C------------------------------------------
1001C NITER_CHL = 0
1002C IF (LAGMOD==1) THEN
1003C ALPHA = LAG_ALPH
1004C IF (LAG_ALPHS==ZERO) THEN
1005C DO IC=1,NC
1006C DIAG_L(IC) = DIAG_H(IC)
1007C ENDDO
1008C ELSE
1009C DO IC=1,NC
1010C DIAG_L(IC) = DIAG_H(IC) + LAG_ALPHS
1011C ENDDO
1012C ENDIF
1013C IP = 1
1014C DO WHILE (IP>0)
1015C DO IH=1,LENH
1016C HL(IH) = HH(IH)
1017C ENDDO
1018C IP = CHOLFACT(NC,DIAG_L,HL,IADH,JCIH,WORK1,WORK2,WORK3)
1019C NITER_CHL = NITER_CHL+1
1020C ALPHA = ALPHA*2
1021C AS = ALPHA + LAG_ALPHS
1022C IF (AS>ASMAX) THEN
1023C IP = 0
1024C WRITE(IOUT,*) '***WARNING (LAGMULT : FACTORISATION FAILED )'
1025C WRITE(ISTDO,*)'***WARNING (LAGMULT : FACTORISATION FAILED )'
1026C ELSEIF (IP>0) THEN
1027C DO IC=1,NC
1028C DIAG_L(IC) = DIAG_H(IC) + ALPHA + LAG_ALPHS
1029C ENDDO
1030C ENDIF
1031C ENDDO
1032C IF (NITER_CHL>1) print*,'FACTOR ITERATIONS = ',NITER_CHL
1033C ENDIF
1034 END IF ! fin proc0
1035C=======================================================================
1036C MUMPS
1037C=======================================================================
1038C
1039C Init structure
1040C
1041 CALL spmd_mumps_ini(mumps_par,1)
1042 mumps_par%ICNTL(18) = 0 ! donnees centralisees sur p0
1043C flag output
1044 mumps_par%ICNTL(3)=-1 ! ou iout par exemple si sorties voulues
1045 IF(ispmd==0) THEN ! structure de stockage sur p0
1046 mumps_par%N = nc ! valid sur proc0 uniquement
1047 mumps_par%NZ = nc+iadh(nc+1)-1 ! valid sur proc0 uniquement
1048 ALLOCATE(mumps_par%A(mumps_par%NZ))
1049 ALLOCATE(mumps_par%IRN(mumps_par%NZ))
1050 ALLOCATE(mumps_par%JCN(mumps_par%NZ))
1051 DO ic = 1, nc
1052 mumps_par%A(ic) = diag_h(ic)
1053 mumps_par%IRN(ic)=ic
1054 mumps_par%JCN(ic)=ic
1055 END DO
1056 nz = nc
1057 DO ic = 1, nc
1058 DO ih=iadh(ic),iadh(ic+1)-1
1059 nz = nz + 1
1060 mumps_par%A(nz) = hh(ih)
1061 mumps_par%IRN(nz)=ic
1062 mumps_par%JCN(nz)=jcih(ih)
1063C MUMPS_PAR%JCN(NZ)=JDIH(IH)
1064 ENDDO
1065 END DO
1066C
1067C calcul de b (r=b)
1068C
1069 DO ic=1,nc
1070 DO ik=iadll(ic),iadll(ic+1)-1
1071 i = indexlag(lll(ik))
1072 j = jll(ik)
1073 IF (j>3) THEN
1074 j = j-3
1075 r(ic) = r(ic) + xll(ik)*(vr(j,i)/dt12+ar(j,i))
1076 ELSE
1077 r(ic) = r(ic) + xll(ik)*(v(j,i)/dt12+a(j,i))
1078 ENDIF
1079 ENDDO
1080 ENDDO
1081 ALLOCATE(mumps_par%RHS(mumps_par%N))
1082 DO ic=1,nc
1083 mumps_par%RHS(ic) = r(ic)
1084 ENDDO
1085 END IF
1086C
1087C Resolution
1088C
1089 CALL spmd_mumps_exec(mumps_par,1) ! analyse et factorisation
1090C
1091 CALL spmd_mumps_exec(mumps_par,2) ! resolution
1092C
1093C Resultat sur p0
1094C
1095 IF(ispmd==0) THEN
1096 DO ic=1,nc
1097 lambda(ic) = mumps_par%RHS(ic)
1098 ENDDO
1099 END IF
1100C
1101C Liberation de la memoire
1102C
1103 CALL spmd_mumps_deal(mumps_par)
1104CC --- Initialisation du vecteur solution
1105C RNORM = ZERO
1106C DO IC=1,NC
1107C RNORM = RNORM + R(IC)*R(IC)
1108C P(IC) = LAMBDA(IC)
1109C ENDDO
1110CC--- Residu initial R(i) = R(i) - sum(H(i,j)*P(j))
1111C DO IC=1,NC
1112C R(IC) = R(IC) - DIAG_H(IC)*P(IC)
1113C DO IH=IADH(IC),IADH(IC+1)-1
1114C JC = JCIH(IH)
1115C R(IC) = R(IC) - HH(IH)*P(JC)
1116C R(JC) = R(JC) - HH(IH)*P(IC)
1117C ENDDO
1118C ENDDO
1119C=======================================================================
1120C iterations
1121C=======================================================================
1122C ITERMAX = NC
1123C ITER = 0
1124C DO WHILE(ITER<ITERMAX)
1125C ITER = ITER+1
1126C NITER = ITER
1127CC--------------------
1128CC preconditionnement z = inv(H) r
1129CC--------------------
1130C IF (LAGMOD==1) THEN
1131CC--- preconditionnement de Cholesky
1132C CALL PRECHOL(Z,DIAG_L,HL,R,NC,IADH,JCIH)
1133C ELSEIF(LAGMOD==2) THEN
1134CC--- preconditionnement polynomial : H=I-B, z = [I+B] r
1135C DO IC=1,NC
1136C Z(IC) = R(IC) + R(IC)
1137C ENDDO
1138C DO IC=1,NC
1139C Z(IC) = Z(IC) - DIAG_H(IC)*R(IC)
1140C DO IH=IADH(IC),IADH(IC+1)-1
1141C JC = JCIH(IH)
1142C Z(IC) = Z(IC) - HH(IH)*R(JC)
1143C Z(JC) = Z(JC) - HH(IH)*R(IC)
1144C ENDDO
1145C ENDDO
1146C ENDIF
1147CC--------------------
1148C R2 = R2NEW
1149C R2NEW = ZERO
1150C DO IC=1,NC
1151C R2NEW = R2NEW + R(IC)*Z(IC)
1152C Q(IC) = ZERO
1153C ENDDO
1154CC--------------------
1155C IF(ITER==1)THEN
1156C DO IC=1,NC
1157C P(IC) = Z(IC)
1158C ENDDO
1159C ELSE
1160C BETA = R2NEW/R2
1161C DO IC=1,NC
1162C P(IC) = Z(IC) + BETA*P(IC)
1163C ENDDO
1164C ENDIF
1165CC-----------------------------------------------------------------------
1166CC calcul de q = q_ + [H] p
1167CC-----------------------------------------------------------------------
1168C DO IC=1,NC
1169C Q(IC) = Q(IC) + DIAG_H(IC)*P(IC)
1170C DO IH=IADH(IC),IADH(IC+1)-1
1171C JC = JCIH(IH)
1172C Q(IC) = Q(IC) + HH(IH)*P(JC)
1173C Q(JC) = Q(JC) + HH(IH)*P(IC)
1174C ENDDO
1175C ENDDO
1176CC-----------------------------------------------------------------------
1177CC pt [H] p = pt q
1178CC-----------------------------------------------------------------------
1179C PQ = ZERO
1180C DO IC=1,NC
1181C PQ = PQ + P(IC)*Q(IC)
1182C ENDDO
1183CC-----------------------------------------------------------------------
1184C IF(PQ==ZERO)THEN
1185C ITER = ITERMAX
1186C ELSE
1187C ALPHA = R2NEW/PQ
1188C LAG_ERSQ2 = ZERO
1189C DO IC=1,NC
1190C LAMBDA(IC) = LAMBDA(IC) + ALPHA*P(IC)
1191C R(IC) = R(IC) - ALPHA*Q(IC)
1192C LAG_ERSQ2 = LAG_ERSQ2 + R(IC)*R(IC)
1193C ENDDO
1194CC--------- test convergence --------------------------------------------
1195C LAG_ERSQ2 = LAG_ERSQ2/MAX(EM30,RNORM)
1196C IF(LAG_ERSQ2<LAGM_TOL)THEN
1197C ITER = ITERMAX
1198C ENDIF
1199C ENDIF
1200C=======================================================================
1201C fin boucle gradient conjugue
1202C=======================================================================
1203C ENDDO
1204C NITER_GC = NITER
1205C---
1206C a = ao + [M]-1[L]t la
1207C---
1208 IF(ispmd==0) THEN
1209 DO ic=1,ncr
1210 DO ik=iadll(ic),iadll(ic+1)-1
1211 i = indexlag(lll(ik))
1212 j = jll(ik)
1213 xll(ik) = xll(ik)*lambda(ic)
1214 IF(j>3) THEN
1215 j = j-3
1216 ar(j,i) = ar(j,i) - xll(ik)/in(i)
1217 ELSE
1218 a(j,i) = a(j,i) - xll(ik)/mas(i)
1219 ENDIF
1220 ENDDO
1221 ENDDO
1222 END IF
1223#else
1224 write(6,*) "Error: this feature requires the MUMPS library "
1225 call flush(6)
1226 CALL arret(5)
1227#endif
1228C--------------------------------------------
1229 RETURN
cholfact
integer function cholfact(n, diag, h, iadh, jcih, wcol, irow, jcol)
Definition cholfact.F:31
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
spmd_mumps_ini
subroutine spmd_mumps_ini(mumps_par, sym)
Definition imp_spmd.F:498
spmd_mumps_deal
subroutine spmd_mumps_deal(mumps_par)
Definition imp_spmd.F:558
spmd_mumps_exec
subroutine spmd_mumps_exec(mumps_par, itask)
Definition imp_spmd.F:724
lag_mult_hp
subroutine lag_mult_hp(nc, lenh, lhmax, ms, in, diag, hh, iadll, lll, jll, xll, ltsm, iadhf, jcihf, iadh, jcih, rbyl, npbyl, icftag, jcftag, ncf_s, ncf_e, ncr, indexlag)
Definition lag_mult_h.F:320
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339
arret
subroutine arret(nn)
Definition arret.F:87

◆ lag_mult_solv()

subroutine lag_mult_solv ( integer nh,
integer nc,
integer ncr,
a,
v,
mas,
integer, dimension(*) iadll,
integer, dimension(*) lll,
integer, dimension(*) jll,
xll,
integer, dimension(*) iadh,
integer, dimension(*) jcih,
hh,
z,
p,
r,
q,
ltsm,
hl,
diag_h,
diag_l,
work1,
integer, dimension(*) work2,
integer, dimension(*) work3,
lambda,
rbyl,
integer, dimension(nnpby,*) npbyl,
ar,
vr,
in,
integer, dimension(*) iadhf,
integer, dimension(*) jcihf,
integer, dimension(*) icftag,
integer, dimension(*) jcftag,
integer ncf_s,
integer ncf_e )

Definition at line 32 of file lag_mult_solv.F.

41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "param_c.inc"
49#include "lagmult.inc"
50#include "com08_c.inc"
51#include "units_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER NC,NCR,NCF_S,NCF_E,NH
56 INTEGER LLL(*),JLL(*),IADLL(*),IADH(*),JCIH(*),IADHF(*),
57 . JCIHF(*),NPBYL(NNPBY,*),WORK2(*),WORK3(*),
58 . ICFTAG(*),JCFTAG(*)
59C REAL
60 my_real
61 . mas(*),in(*), a(3,*),ar(3,*),v(3,*),vr(3,*),
62 . xll(*),p(*),r(*),q(*),lambda(*),rbyl(nrby,*),
63 . ltsm(6,*),z(*),hh(*),hl(*),diag_h(*),diag_l(*),
64 . work1(*)
65C-----------------------------------------------
66C L o c a l V a r i a b l e s
67C-----------------------------------------------
68 INTEGER I,J,IK,IC,JC,IH,IP,ITER,NITER,ITERMAX,NIN,LENH
69 my_real
70 . as,asmax,alpha,beta,r2,r2new,rnorm,pq,hij,scale,dd
71C-----------------------------------------------
72C External Functions
73C-----------------------------------------------
74 INTEGER CHOLFACT
75 EXTERNAL cholfact
76C-----------------------------------------------
77C
78C | M | Lt| | a | M ao
79C |---+---| | = |
80C | L | 0 | | la | bo
81C
82C [M] a + [L]t la = [M] ao
83C [L] a = bo
84C
85C a = -[M]-1[L]t la + ao
86C [L][M]-1[L]t la = [L] ao - bo
87C
88C [H] = [L][M]-1[L]t
89C b = [L] ao - bo
90C [H] la = b
91C
92C a = ao - [M]-1[L]t la
93C-----------------------------------------------
94C
95C la : LAMBDA(NC)
96C ao : A(3,NUMNOD)
97C [M](i) : MAS(NUMNOD)
98C [L][M]-1[L]t la : Q(NC)
99C [L] ao - b : B(NC)
100C [M]-1[L]t : LTSM(3,NUMNOD) une colonne de [M]-1[L]t
101C : Z(NC)
102C
103C [L](ic,j,i) -> XLL(IK)
104C ic = 1..nc : IK = IADLL(IC)...IADLL(IC+1)-1
105C j = 1,2,3 : J = JLL(IK)
106C i = 1..numnod : I = LLL(IK)
107C
108C
109C NC : nombre de contact
110C
111C IC : numero du contact (1,NC)
112C IK : index de XLL,LLL,JLL
113C I : numero global du noeud (1,NUMNOD)
114C J : direction 1,2,3
115C
116C-----------------------------------------------
117C [H](ic,jc) -> HH(IH)
118C ic = 1..nc : IH = IADH(IC)...IADH(IC+1)-1
119C jc = 1..ic : JC = JCIH(IH)
120C
121C q = [H] p
122C
123C do ic=1,nc
124C q(ic) = 0
125C do jc=1,nc
126C q(ic) = q(ic) + h(ic,jc)*p(jc)
127C enddo
128C enddo
129C
130C DO IC=1,NC
131C Q(IC) = 0
132C DO IH=IADH(IC),IADH(IC+1)-1
133C Q(IC) = Q(IC) - HH(IH)*P(JCIH(IH))
134C ENDDO
135C ENDDO
136C-----------------------------------------------
137C evaluation de b:
138C
139C Vs = Somme(Ni Vi) + vout
140C Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai) + vout
141C Somme(dt Ni Ai) - dt As = Vs_ -Somme(Ni Vi_) - vout
142C [L] = dt {N1,N2,..,N15,-1}
143C bo = {N1,N2,..,N15,-1} a = -[L]/dt v_ - vout
144C b = [L] ao - bo
145C b = [L] ao + [L]/dt v_ + vout
146C b = [L]/dt vo+ + vout
147C-----------------------------------------------------------------------
148C
149C calcul de [H] la = b par gradient conjugue
150C
151C-----------------------------------------------------------------------
152C Gradient conjugue [H] la = b avec preconditionnement [N]
153C [H]-1 = Somme([I-H]^i),i=0,1...inf si |I-H|<1
154C [N]-1 = [2I-H]
155C-----------------------------------------------------------------------
156C r(0) = b - [H] la(0)
157C for i=1,niter
158C z(i-1) = [N]-1 r(i-1)
159C rho(i-1) = r(i-1)T z(i-1)
160C if (i == 1)
161C p(1) = z(0)
162C else
163C beta = rho(i-1) / rho(i-2)
164C p(i) = z(i-1) + beta p(i-1)
165C endif
166C q(i) = [H] p(i)
167C alpha = rho(i-1) / p(i)T q(i)
168C la(i) = la(i-1) + alpha p(i)
169C r(i) = r(i-1) - alpha q(i)
170C end
171C======================================================================|
172 scale = three_over_4
173 asmax = zep6
174 r2new = zero
175 lag_nc= nc
176C---
177C--- Calcul de la matrice H
178 CALL lag_mult_h(nc ,lenh ,nh ,mas ,in ,
179 2 diag_h ,hh ,iadll ,lll ,jll ,
180 3 xll ,ltsm ,iadhf ,jcihf ,iadh ,
181 4 jcih ,rbyl ,npbyl ,icftag ,jcftag ,
182 5 ncf_s ,ncf_e ,ncr )
183C------------------------------------------
184C Normalisation du systeme
185c lagopt=0 => sans norme
186c lagopt=1 => norme L1 (diagonal)
187c lagopt=2 => norme L2 (columns)
188C------------------------------------------
189 IF (lagopt==1) THEN
190 dd = zero
191 DO ic=1,nc
192 dd = dd + diag_h(ic)*diag_h(ic)
193 ENDDO
194 DO ic=1,nc
195 z(ic) = dd
196 ENDDO
197 ELSEIF (lagopt==2) THEN
198 DO ic=1,nc
199 z(ic) = diag_h(ic)*diag_h(ic)
200 ENDDO
201 DO ic=1,nc
202 DO ih=iadh(ic),iadh(ic+1)-1
203 jc = jcih(ih)
204 hij = hh(ih)*hh(ih)
205 z(ic) = z(ic) + hij
206 z(jc) = z(jc) + hij
207 ENDDO
208 ENDDO
209 ENDIF
210C---
211 IF (lagopt>0) THEN
212 DO ic=1,nc
213 z(ic) = scale/sqrt(sqrt(z(ic)))
214 ENDDO
215C--- normalisation [H] Somme([H](ic,i)^2)=1 i=1...nc
216 DO ic=1,nc
217 diag_h(ic) = diag_h(ic)*z(ic)*z(ic)
218 DO ih=iadh(ic),iadh(ic+1)-1
219 hh(ih) = hh(ih)*z(ic)*z(jcih(ih))
220 ENDDO
221 ENDDO
222C--- normalisation [LLL] et {R}
223 DO ic=1,nc
224 DO ik=iadll(ic),iadll(ic+1)-1
225 xll(ik) = xll(ik)*z(ic)
226 ENDDO
227 r(ic) = r(ic)*z(ic)
228 ENDDO
229 ENDIF
230C------------------------------------------
231C Factorisation de [H] : Cholesky incomplet + shift iteratif
232C------------------------------------------
233 niter_chl = 0
234 IF (lagmod==1) THEN
235 alpha = lag_alph
236 IF (lag_alphs==zero) THEN
237 DO ic=1,nc
238 diag_l(ic) = diag_h(ic)
239 ENDDO
240 ELSE
241 DO ic=1,nc
242 diag_l(ic) = diag_h(ic) + lag_alphs
243 ENDDO
244 ENDIF
245 ip = 1
246 DO WHILE (ip>0)
247 DO ih=1,lenh
248 hl(ih) = hh(ih)
249 ENDDO
250 ip = cholfact(nc,diag_l,hl,iadh,jcih,work1,work2,work3)
251 niter_chl = niter_chl+1
252 alpha = alpha*2
253 as = alpha + lag_alphs
254 IF (as>asmax) THEN
255 ip = 0
256 WRITE(iout,*) '***WARNING (LAGMULT : FACTORISATION FAILED )'
257 WRITE(istdo,*)'***WARNING (LAGMULT : FACTORISATION FAILED )'
258 ELSEIF (ip>0) THEN
259 DO ic=1,nc
260 diag_l(ic) = diag_h(ic) + alpha + lag_alphs
261 ENDDO
262 ENDIF
263 ENDDO
264 IF (niter_chl>1) print*,'FACTOR ITERATIONS = ',niter_chl
265 ENDIF
266C=======================================================================
267C calcul de b (r=b)
268C=======================================================================
269 DO ic=1,nc
270 DO ik=iadll(ic),iadll(ic+1)-1
271 i = lll(ik)
272 j = jll(ik)
273 IF (j>3) THEN
274 j = j-3
275 r(ic) = r(ic) + xll(ik)*(vr(j,i)/dt12+ar(j,i))
276 ELSE
277 r(ic) = r(ic) + xll(ik)*(v(j,i)/dt12+a(j,i))
278 ENDIF
279 ENDDO
280 ENDDO
281C--- Initialisation du vecteur solution
282 rnorm = zero
283 DO ic=1,nc
284 rnorm = rnorm + r(ic)*r(ic)
285 p(ic) = lambda(ic)
286 ENDDO
287C--- Residu initial R(i) = R(i) - sum(H(i,j)*P(j))
288 DO ic=1,nc
289 r(ic) = r(ic) - diag_h(ic)*p(ic)
290 DO ih=iadh(ic),iadh(ic+1)-1
291 jc = jcih(ih)
292 r(ic) = r(ic) - hh(ih)*p(jc)
293 r(jc) = r(jc) - hh(ih)*p(ic)
294 ENDDO
295 ENDDO
296C=======================================================================
297C iterations
298C=======================================================================
299 itermax = nc
300 iter = 0
301 niter = 0
302 DO WHILE(iter<itermax)
303 iter = iter+1
304 niter = iter
305C--------------------
306C preconditionnement z = inv(H) r
307C--------------------
308 IF (lagmod==1) THEN
309C--- preconditionnement de Cholesky
310 CALL prechol(z,diag_l,hl,r,nc,iadh,jcih)
311 ELSEIF(lagmod==2) THEN
312C--- preconditionnement polynomial : H=I-B, z = [I+B] r
313 DO ic=1,nc
314 z(ic) = r(ic) + r(ic)
315 ENDDO
316 DO ic=1,nc
317 z(ic) = z(ic) - diag_h(ic)*r(ic)
318 DO ih=iadh(ic),iadh(ic+1)-1
319 jc = jcih(ih)
320 z(ic) = z(ic) - hh(ih)*r(jc)
321 z(jc) = z(jc) - hh(ih)*r(ic)
322 ENDDO
323 ENDDO
324 ENDIF
325C--------------------
326 r2 = r2new
327 r2new = zero
328 DO ic=1,nc
329 r2new = r2new + r(ic)*z(ic)
330 q(ic) = zero
331 ENDDO
332C--------------------
333 IF(iter==1)THEN
334 DO ic=1,nc
335 p(ic) = z(ic)
336 ENDDO
337 ELSE
338 beta = r2new/r2
339 DO ic=1,nc
340 p(ic) = z(ic) + beta*p(ic)
341 ENDDO
342 ENDIF
343C-----------------------------------------------------------------------
344C calcul de q = q_ + [H] p
345C-----------------------------------------------------------------------
346 DO ic=1,nc
347 q(ic) = q(ic) + diag_h(ic)*p(ic)
348 DO ih=iadh(ic),iadh(ic+1)-1
349 jc = jcih(ih)
350 q(ic) = q(ic) + hh(ih)*p(jc)
351 q(jc) = q(jc) + hh(ih)*p(ic)
352 ENDDO
353 ENDDO
354C-----------------------------------------------------------------------
355C pt [H] p = pt q
356C-----------------------------------------------------------------------
357 pq = zero
358 DO ic=1,nc
359 pq = pq + p(ic)*q(ic)
360 ENDDO
361C-----------------------------------------------------------------------
362 IF(pq==zero)THEN
363 iter = itermax
364 ELSE
365 alpha = r2new/pq
366 lag_ersq2 = zero
367 DO ic=1,nc
368 lambda(ic) = lambda(ic) + alpha*p(ic)
369 r(ic) = r(ic) - alpha*q(ic)
370 lag_ersq2 = lag_ersq2 + r(ic)*r(ic)
371 ENDDO
372C--------- test convergence --------------------------------------------
373 lag_ersq2 = lag_ersq2/max(em30,rnorm)
374 IF(lag_ersq2<lagm_tol)THEN
375 iter = itermax
376 ENDIF
377 ENDIF
378C=======================================================================
379C fin boucle gradient conjugue
380C=======================================================================
381 ENDDO
382 niter_gc = niter
383C---
384C a = ao + [M]-1[L]t la
385C---
386 DO ic=1,ncr
387 DO ik=iadll(ic),iadll(ic+1)-1
388 i = lll(ik)
389 j = jll(ik)
390 xll(ik) = xll(ik)*lambda(ic)
391 IF(j>3) THEN
392 j = j-3
393 ar(j,i) = ar(j,i) - xll(ik)/in(i)
394 ELSE
395 a(j,i) = a(j,i) - xll(ik)/mas(i)
396 ENDIF
397 ENDDO
398 ENDDO
399C--------------------------------------------
400 RETURN
prechol
subroutine prechol(z, d, l, r, nc, iadh, jcih)
Definition cholfact.F:119
lag_mult_h
subroutine lag_mult_h(nc, lenh, lhmax, ms, in, diag, hh, iadll, lll, jll, xll, ltsm, iadhf, jcihf, iadh, jcih, rbyl, npbyl, icftag, jcftag, ncf_s, ncf_e, ncr)
Definition lag_mult_h.F:39
max
#define max(a, b)
Definition macros.h:21

◆ lag_mult_solvp()

subroutine lag_mult_solvp ( integer nh,
integer nc,
integer ncr,
a,
v,
mas,
integer, dimension(*) iadll,
integer, dimension(*) lll,
integer, dimension(*) jll,
xll,
integer, dimension(*) iadh,
integer, dimension(*) jcih,
hh,
z,
p,
r,
q,
ltsm,
hl,
diag_h,
diag_l,
work1,
integer, dimension(*) work2,
integer, dimension(*) work3,
lambda,
rbyl,
integer, dimension(nnpby,*) npbyl,
ar,
vr,
in,
integer, dimension(*) iadhf,
integer, dimension(*) jcihf,
integer, dimension(*) icftag,
integer, dimension(*) jcftag,
integer ncf_s,
integer ncf_e,
integer, dimension(*) indexlag )

Definition at line 412 of file lag_mult_solv.F.

421C-----------------------------------------------
422C I m p l i c i t T y p e s
423C-----------------------------------------------
424#include "implicit_f.inc"
425C-----------------------------------------------
426C C o m m o n B l o c k s
427C-----------------------------------------------
428#include "param_c.inc"
429#include "lagmult.inc"
430#include "com08_c.inc"
431#include "units_c.inc"
432C-----------------------------------------------
433C D u m m y A r g u m e n t s
434C-----------------------------------------------
435 INTEGER NC,NCR,NCF_S,NCF_E,NH
436 INTEGER LLL(*),JLL(*),IADLL(*),IADH(*),JCIH(*),IADHF(*),
437 . JCIHF(*),NPBYL(NNPBY,*),WORK2(*),WORK3(*),
438 . ICFTAG(*),JCFTAG(*), INDEXLAG(*)
439C REAL
440 my_real
441 . mas(*),in(*), a(3,*),ar(3,*),v(3,*),vr(3,*),
442 . xll(*),p(*),r(*),q(*),lambda(*),rbyl(nrby,*),
443 . ltsm(6,*),z(*),hh(*),hl(*),diag_h(*),diag_l(*),
444 . work1(*)
445C-----------------------------------------------
446C L o c a l V a r i a b l e s
447C-----------------------------------------------
448 INTEGER I,J,IK,IC,JC,IH,IP,ITER,NITER,ITERMAX,NIN,LENH,
449 . NZ
450 my_real
451 . as,asmax,alpha,beta,r2,r2new,rnorm,pq,hij,scale,dd
452C-----------------------------------------------
453C External Functions
454C-----------------------------------------------
455 INTEGER CHOLFACT
456 EXTERNAL cholfact
457C-----------------------------------------------
458C
459C | M | Lt| | a | M ao
460C |---+---| | = |
461C | L | 0 | | la | bo
462C
463C [M] a + [L]t la = [M] ao
464C [L] a = bo
465C
466C a = -[M]-1[L]t la + ao
467C [L][M]-1[L]t la = [L] ao - bo
468C
469C [H] = [L][M]-1[L]t
470C b = [L] ao - bo
471C [H] la = b
472C
473C a = ao - [M]-1[L]t la
474C-----------------------------------------------
475C
476C la : LAMBDA(NC)
477C ao : A(3,NUMNOD)
478C [M](i) : MAS(NUMNOD)
479C [L][M]-1[L]t la : Q(NC)
480C [L] ao - b : B(NC)
481C [M]-1[L]t : LTSM(3,NUMNOD) une colonne de [M]-1[L]t
482C : Z(NC)
483C
484C [L](ic,j,i) -> XLL(IK)
485C ic = 1..nc : IK = IADLL(IC)...IADLL(IC+1)-1
486C j = 1,2,3 : J = JLL(IK)
487C i = 1..numnod : I = LLL(IK)
488C
489C
490C NC : nombre de contact
491C
492C IC : numero du contact (1,NC)
493C IK : index de XLL,LLL,JLL
494C I : numero global du noeud (1,NUMNOD)
495C J : direction 1,2,3
496C
497C-----------------------------------------------
498C [H](ic,jc) -> HH(IH)
499C ic = 1..nc : IH = IADH(IC)...IADH(IC+1)-1
500C jc = 1..ic : JC = JCIH(IH)
501C
502C q = [H] p
503C
504C do ic=1,nc
505C q(ic) = 0
506C do jc=1,nc
507C q(ic) = q(ic) + h(ic,jc)*p(jc)
508C enddo
509C enddo
510C
511C DO IC=1,NC
512C Q(IC) = 0
513C DO IH=IADH(IC),IADH(IC+1)-1
514C Q(IC) = Q(IC) - HH(IH)*P(JCIH(IH))
515C ENDDO
516C ENDDO
517C-----------------------------------------------
518C evaluation de b:
519C
520C Vs = Somme(Ni Vi) + vout
521C Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai) + vout
522C Somme(dt Ni Ai) - dt As = Vs_ -Somme(Ni Vi_) - vout
523C [L] = dt {N1,N2,..,N15,-1}
524C bo = {N1,N2,..,N15,-1} a = -[L]/dt v_ - vout
525C b = [L] ao - bo
526C b = [L] ao + [L]/dt v_ + vout
527C b = [L]/dt vo+ + vout
528C-----------------------------------------------------------------------
529C
530C calcul de [H] la = b par gradient conjugue
531C
532C-----------------------------------------------------------------------
533C Gradient conjugue [H] la = b avec preconditionnement [N]
534C [H]-1 = Somme([I-H]^i),i=0,1...inf si |I-H|<1
535C [N]-1 = [2I-H]
536C-----------------------------------------------------------------------
537C r(0) = b - [H] la(0)
538C for i=1,niter
539C z(i-1) = [N]-1 r(i-1)
540C rho(i-1) = r(i-1)T z(i-1)
541C if (i == 1)
542C p(1) = z(0)
543C else
544C beta = rho(i-1) / rho(i-2)
545C p(i) = z(i-1) + beta p(i-1)
546C endif
547C q(i) = [H] p(i)
548C alpha = rho(i-1) / p(i)T q(i)
549C la(i) = la(i-1) + alpha p(i)
550C r(i) = r(i-1) - alpha q(i)
551C end
552C======================================================================|
553 scale = three_over_4
554 asmax = zep6
555 r2new = zero
556 lag_nc= nc
557C---
558C--- Calcul de la matrice H
559 CALL lag_mult_hp(nc ,lenh ,nh ,mas ,in ,
560 2 diag_h ,hh ,iadll ,lll ,jll ,
561 3 xll ,ltsm ,iadhf ,jcihf ,iadh ,
562 4 jcih ,rbyl ,npbyl ,icftag ,jcftag ,
563 5 ncf_s ,ncf_e ,ncr ,indexlag)
564C------------------------------------------
565C Normalisation du systeme
566c lagopt=0 => sans norme
567c lagopt=1 => norme L1 (diagonal)
568c lagopt=2 => norme L2 (columns)
569C------------------------------------------
570 IF (lagopt==1) THEN
571 dd = zero
572 DO ic=1,nc
573 dd = dd + diag_h(ic)*diag_h(ic)
574 ENDDO
575 DO ic=1,nc
576 z(ic) = dd
577 ENDDO
578 ELSEIF (lagopt==2) THEN
579 DO ic=1,nc
580 z(ic) = diag_h(ic)*diag_h(ic)
581 ENDDO
582 DO ic=1,nc
583 DO ih=iadh(ic),iadh(ic+1)-1
584 jc = jcih(ih)
585 hij = hh(ih)*hh(ih)
586 z(ic) = z(ic) + hij
587 z(jc) = z(jc) + hij
588 ENDDO
589 ENDDO
590 ENDIF
591C---
592 IF (lagopt>0) THEN
593 DO ic=1,nc
594 z(ic) = scale/sqrt(sqrt(z(ic)))
595 ENDDO
596C--- normalisation [H] Somme([H](ic,i)^2)=1 i=1...nc
597 DO ic=1,nc
598 diag_h(ic) = diag_h(ic)*z(ic)*z(ic)
599 DO ih=iadh(ic),iadh(ic+1)-1
600 hh(ih) = hh(ih)*z(ic)*z(jcih(ih))
601 ENDDO
602 ENDDO
603C--- normalisation [LLL] et {R}
604 DO ic=1,nc
605 DO ik=iadll(ic),iadll(ic+1)-1
606 xll(ik) = xll(ik)*z(ic)
607 ENDDO
608 r(ic) = r(ic)*z(ic)
609 ENDDO
610 ENDIF
611C------------------------------------------
612C Factorisation de [H] : Cholesky incomplet + shift iteratif
613C------------------------------------------
614 niter_chl = 0
615 IF (lagmod==1) THEN
616 alpha = lag_alph
617 IF (lag_alphs==zero) THEN
618 DO ic=1,nc
619 diag_l(ic) = diag_h(ic)
620 ENDDO
621 ELSE
622 DO ic=1,nc
623 diag_l(ic) = diag_h(ic) + lag_alphs
624 ENDDO
625 ENDIF
626 ip = 1
627 DO WHILE (ip>0)
628 DO ih=1,lenh
629 hl(ih) = hh(ih)
630 ENDDO
631 ip = cholfact(nc,diag_l,hl,iadh,jcih,work1,work2,work3)
632 niter_chl = niter_chl+1
633 alpha = alpha*2
634 as = alpha + lag_alphs
635 IF (as>asmax) THEN
636 ip = 0
637 WRITE(iout,*) '***WARNING (LAGMULT : FACTORISATION FAILED )'
638 WRITE(istdo,*)'***WARNING (LAGMULT : FACTORISATION FAILED )'
639 ELSEIF (ip>0) THEN
640 DO ic=1,nc
641 diag_l(ic) = diag_h(ic) + alpha + lag_alphs
642 ENDDO
643 ENDIF
644 ENDDO
645 IF (niter_chl>1) print*,'FACTOR ITERATIONS = ',niter_chl
646 ENDIF
647C=======================================================================
648C calcul de b (r=b)
649C=======================================================================
650 DO ic=1,nc
651 DO ik=iadll(ic),iadll(ic+1)-1
652 i = indexlag(lll(ik))
653 j = jll(ik)
654 IF (j>3) THEN
655 j = j-3
656 r(ic) = r(ic) + xll(ik)*(vr(j,i)/dt12+ar(j,i))
657 ELSE
658 r(ic) = r(ic) + xll(ik)*(v(j,i)/dt12+a(j,i))
659 ENDIF
660 ENDDO
661 ENDDO
662C--- Initialisation du vecteur solution
663 rnorm = zero
664 DO ic=1,nc
665 rnorm = rnorm + r(ic)*r(ic)
666 p(ic) = lambda(ic)
667 ENDDO
668C--- Residu initial R(i) = R(i) - sum(H(i,j)*P(j))
669 DO ic=1,nc
670 r(ic) = r(ic) - diag_h(ic)*p(ic)
671 DO ih=iadh(ic),iadh(ic+1)-1
672 jc = jcih(ih)
673 r(ic) = r(ic) - hh(ih)*p(jc)
674 r(jc) = r(jc) - hh(ih)*p(ic)
675 ENDDO
676 ENDDO
677C=======================================================================
678C iterations
679C=======================================================================
680 itermax = nc
681 iter = 0
682 DO WHILE(iter<itermax)
683 iter = iter+1
684 niter = iter
685C--------------------
686C preconditionnement z = inv(H) r
687C--------------------
688 IF (lagmod==1) THEN
689C--- preconditionnement de Cholesky
690 CALL prechol(z,diag_l,hl,r,nc,iadh,jcih)
691 ELSEIF(lagmod==2) THEN
692C--- preconditionnement polynomial : H=I-B, z = [I+B] r
693 DO ic=1,nc
694 z(ic) = r(ic) + r(ic)
695 ENDDO
696 DO ic=1,nc
697 z(ic) = z(ic) - diag_h(ic)*r(ic)
698 DO ih=iadh(ic),iadh(ic+1)-1
699 jc = jcih(ih)
700 z(ic) = z(ic) - hh(ih)*r(jc)
701 z(jc) = z(jc) - hh(ih)*r(ic)
702 ENDDO
703 ENDDO
704 ENDIF
705C--------------------
706 r2 = r2new
707 r2new = zero
708 DO ic=1,nc
709 r2new = r2new + r(ic)*z(ic)
710 q(ic) = zero
711 ENDDO
712C--------------------
713 IF(iter==1)THEN
714 DO ic=1,nc
715 p(ic) = z(ic)
716 ENDDO
717 ELSE
718 beta = r2new/r2
719 DO ic=1,nc
720 p(ic) = z(ic) + beta*p(ic)
721 ENDDO
722 ENDIF
723C-----------------------------------------------------------------------
724C calcul de q = q_ + [H] p
725C-----------------------------------------------------------------------
726 DO ic=1,nc
727 q(ic) = q(ic) + diag_h(ic)*p(ic)
728 DO ih=iadh(ic),iadh(ic+1)-1
729 jc = jcih(ih)
730 q(ic) = q(ic) + hh(ih)*p(jc)
731 q(jc) = q(jc) + hh(ih)*p(ic)
732 ENDDO
733 ENDDO
734C-----------------------------------------------------------------------
735C pt [H] p = pt q
736C-----------------------------------------------------------------------
737 pq = zero
738 DO ic=1,nc
739 pq = pq + p(ic)*q(ic)
740 ENDDO
741C-----------------------------------------------------------------------
742 IF(pq==zero)THEN
743 iter = itermax
744 ELSE
745 alpha = r2new/pq
746 lag_ersq2 = zero
747 DO ic=1,nc
748 lambda(ic) = lambda(ic) + alpha*p(ic)
749 r(ic) = r(ic) - alpha*q(ic)
750 lag_ersq2 = lag_ersq2 + r(ic)*r(ic)
751 ENDDO
752C--------- test convergence --------------------------------------------
753 lag_ersq2 = lag_ersq2/max(em30,rnorm)
754 IF(lag_ersq2<lagm_tol)THEN
755 iter = itermax
756 ENDIF
757 ENDIF
758C=======================================================================
759C fin boucle gradient conjugue
760C=======================================================================
761 ENDDO
762 niter_gc = niter
763C---
764C a = ao + [M]-1[L]t la
765C---
766 DO ic=1,ncr
767 DO ik=iadll(ic),iadll(ic+1)-1
768 i = indexlag(lll(ik))
769 j = jll(ik)
770 xll(ik) = xll(ik)*lambda(ic)
771 IF(j>3) THEN
772 j = j-3
773 ar(j,i) = ar(j,i) - xll(ik)/in(i)
774 ELSE
775 a(j,i) = a(j,i) - xll(ik)/mas(i)
776 ENDIF
777 ENDDO
778 ENDDO
779C--------------------------------------------
780 RETURN