prebd4.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\prebd4.f
00002 !
00059                      SUBROUTINE PREBD4
00060 !                    *****************
00061 !
00062      &(X1,X2,A11,A12,A21,A22,B1,B2,D11,D12,D21,D22,
00063      & MESH,PREXSM,DIADON)
00064 !
00065 !***********************************************************************
00066 ! BIEF   V6P1                                   21/08/2010
00067 !***********************************************************************
00068 !
00069 !
00070 !
00071 !
00072 !
00073 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00074 !| A11            |<->| TERM (1,1) OF MATRIX
00075 !| A12            |<->| TERM (1,2) OF MATRIX
00076 !| A21            |<->| TERM (2,1) OF MATRIX
00077 !| A22            |<->| TERM (2,2) OF MATRIX
00078 !| B1             |<->| FIRST RIGHT-HAND SIDE
00079 !| B2             |<->| SECOND RIGHT-HAND SIDE
00080 !| D11            |<--| DIAGONAL MATRIX
00081 !| D12            |<--| DIAGONAL MATRIX
00082 !| D21            |<--| DIAGONAL MATRIX
00083 !| D22            |<--| DIAGONAL MATRIX
00084 !| DIADON         |-->| .TRUE. : DIAGONALS ARE GIVEN
00085 !| MESH           |-->| MESH STRUCTURE
00086 !| PREXSM         |-->| .TRUE. : PRECONDITIONING X1,X2 AND B1,B2
00087 !| X1             |<->| FIRST INITIAL GUESS
00088 !| X2             |-->| SECOND INITIAL GUESS
00089 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00090 !
00091       USE BIEF, EX_PREBD4 => PREBD4
00092 !
00093       IMPLICIT NONE
00094       INTEGER LNG,LU
00095       COMMON/INFO/LNG,LU
00096 !
00097 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00098 !
00099       LOGICAL, INTENT(IN) :: PREXSM,DIADON
00100 !
00101 !-----------------------------------------------------------------------
00102 !
00103 !  VECTOR STRUCTURES
00104 !
00105       TYPE(BIEF_OBJ), INTENT(INOUT) :: X1,B2,D11,D12,D21,D22
00106       TYPE(BIEF_OBJ), INTENT(IN)    :: X2,B1
00107 !
00108 !-----------------------------------------------------------------------
00109 !
00110 !  MATRIX STRUCTURES
00111 !
00112       TYPE(BIEF_OBJ), INTENT(INOUT) :: A11,A12,A21,A22
00113 !
00114 !-----------------------------------------------------------------------
00115 !
00116 !  MESH STRUCTURE
00117 !
00118       TYPE(BIEF_MESH), INTENT(INOUT) :: MESH
00119 !
00120 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00121 !
00122       INTEGER I,NPOIN1,NPOIN2
00123       DOUBLE PRECISION C
00124 !
00125 !-----------------------------------------------------------------------
00126 !
00127       NPOIN1 = X1%DIM1
00128       NPOIN2 = X2%DIM1
00129 !
00130       IF(NPOIN2.NE.NPOIN1) THEN
00131         IF(LNG.EQ.1) WRITE(LU,100)
00132         IF(LNG.EQ.2) WRITE(LU,200)
00133 100     FORMAT(1X,'PREBD4 (BIEF) : MATRICES RECTANGULAIRES',/,1X,
00134      &  'PRECONDITIONNEMENT BLOC-DIAGONAL IMPOSSIBLE DANS CE CAS')
00135 200     FORMAT(1X,'PREBD4 (BIEF) : RECTANGULAR MATRICES',/,1X,
00136      &  'BLOCK-DIAGONAL PRECONDITIONING IMPOSSIBLE IN THIS CASE')
00137         CALL PLANTE(1)
00138         STOP
00139       ENDIF
00140 !
00141 !-----------------------------------------------------------------------
00142 !
00143 !  PREPARES THE DIAGONALS:
00144 !
00145       IF(.NOT.DIADON) THEN
00146 !
00147         CALL OS( 'X=Y     ' , D11 , A11%D , D11 , C )
00148         CALL OS( 'X=Y     ' , D12 , A12%D , D12 , C )
00149         CALL OS( 'X=Y     ' , D21 , A21%D , D21 , C )
00150         CALL OS( 'X=Y     ' , D22 , A22%D , D22 , C )
00151 !
00152 !  TEST TO REDUCE TO DIAGONAL PRECONDITIONING
00153 !
00154 !       CALL OS( 'X=C     ' , D12 , A12%D , Z , 0.D0 )
00155 !       CALL OS( 'X=C     ' , D21 , A21%D , Z , 0.D0 )
00156 !
00157 !  END OF TEST TO REDUCE TO DIAGONAL PRECONDITIONING
00158 !
00159       ENDIF
00160 !
00161 !-----------------------------------------------------------------------
00162 !
00163 !  L D U FACTORISATION OF THE DIAGONAL BLOCK:
00164 !
00165 !     ONLY D11 INVERTED IS NOW USED
00166       CALL OS( 'X=1/Y   ' , D11 , D11 , D11 , C)
00167 !
00168       DO I = 1,NPOIN1
00169 !
00170         D21%R(I) = D21%R(I) * D11%R(I)
00171         D22%R(I) = D22%R(I) - D21%R(I) * D12%R(I)
00172         D12%R(I) = D12%R(I) * D11%R(I)
00173 !
00174       ENDDO
00175 !
00176 !-----------------------------------------------------------------------
00177 !
00178 ! CHANGE OF VARIABLES:
00179 !
00180       IF(PREXSM) THEN
00181 !
00182         CALL OS( 'X=X+YZ  ' , X1 , X2 , D12 , C )
00183 !
00184       ENDIF
00185 !
00186 !  COMPUTES THE SQUARE ROOT
00187 !  INVERTS D11,D22,D33
00188 !  (THEY ARE ONLY USED IN THIS FORM FROM NOW ON)
00189 !
00190 !     INVERSION OF D11 ALREADY PERFORMED
00191       CALL OS( 'X=1/Y   ' , D22 , D22 , D22 , C )
00192       CALL OS( 'X=SQR(Y)' , D11 , D11 , D11 , C )
00193       CALL OS( 'X=SQR(Y)' , D22 , D22 , D22 , C )
00194 !
00195 !=======================================================================
00196 ! MULTIPLIES A ON THE LEFT BY L INVERTED
00197 !=======================================================================
00198 !
00199 ! A21 :
00200       CALL OM( 'M=M-DN  ' , A21 , A11 , D21 , C , MESH)
00201 ! A22 :
00202       CALL OM( 'M=M-DN  ' , A22 , A12 , D21 , C , MESH)
00203 !
00204 !=======================================================================
00205 ! MULTIPLIES A ON THE RIGHT BY U INVERTED
00206 !=======================================================================
00207 !
00208 ! A12 :
00209       CALL OM( 'M=M-ND  ' , A12 , A11 , D12 , C , MESH)
00210 ! A22 :
00211       CALL OM( 'M=M-ND  ' , A22 , A21 , D12 , C , MESH)
00212 !
00213 !-----------------------------------------------------------------------
00214 !
00215 ! NEW SECOND MEMBER
00216 !
00217       IF(PREXSM) THEN
00218 !
00219       DO I = 1,NPOIN1
00220         B2%R(I) = B2%R(I) - D21%R(I) * B1%R(I)
00221       ENDDO
00222 !
00223       ENDIF
00224 !
00225 !-----------------------------------------------------------------------
00226 !
00227       RETURN
00228       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0