diff --git a/src/core_init_atmosphere/mpas_atm_advection.F b/src/core_init_atmosphere/mpas_atm_advection.F index 7d9c7beb24..f2cb3b3064 100644 --- a/src/core_init_atmosphere/mpas_atm_advection.F +++ b/src/core_init_atmosphere/mpas_atm_advection.F @@ -13,6 +13,7 @@ module atm_advection use mpas_constants use mpas_abort, only : mpas_dmpar_global_abort use mpas_log, only : mpas_log_write + use mpas_matrix_operations, only : mpas_migs contains @@ -602,7 +603,7 @@ subroutine poly_fit_2(a_in,b_out,weights_in,m,n,ne) ! call migs(a,n,b,indx) ! else - call migs(atha,n,atha_inv,indx) + call mpas_migs(atha,n,atha_inv,indx) b = matmul(atha_inv,ath) @@ -612,135 +613,7 @@ subroutine poly_fit_2(a_in,b_out,weights_in,m,n,ne) b_out(1:n,1:m) = b(1:n,1:m) end subroutine poly_fit_2 - - - ! Updated 10/24/2001. - ! - !!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ! - ! Please Note: ! - ! ! - ! (1) This computer program is written by Tao Pang in conjunction with ! - ! his book, "An Introduction to Computational Physics," published ! - ! by Cambridge University Press in 1997. ! - ! ! - ! (2) No warranties, express or implied, are made for this program. ! - ! ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - SUBROUTINE MIGS (A,N,X,INDX) - ! - ! Subroutine to invert matrix A(N,N) with the inverse stored - ! in X(N,N) in the output. Copyright (c) Tao Pang 2001. - ! - IMPLICIT NONE - INTEGER, INTENT (IN) :: N - INTEGER :: I,J,K - INTEGER, INTENT (OUT), DIMENSION (N) :: INDX - REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A - REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X - REAL (kind=RKIND), DIMENSION (N,N) :: B - ! - DO I = 1, N - DO J = 1, N - B(I,J) = 0.0 - END DO - END DO - DO I = 1, N - B(I,I) = 1.0 - END DO - ! - CALL ELGS (A,N,INDX) - ! - DO I = 1, N-1 - DO J = I+1, N - DO K = 1, N - B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K) - END DO - END DO - END DO - ! - DO I = 1, N - X(N,I) = B(INDX(N),I)/A(INDX(N),N) - DO J = N-1, 1, -1 - X(J,I) = B(INDX(J),I) - DO K = J+1, N - X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I) - END DO - X(J,I) = X(J,I)/A(INDX(J),J) - END DO - END DO - END SUBROUTINE MIGS - - - SUBROUTINE ELGS (A,N,INDX) - ! - ! Subroutine to perform the partial-pivoting Gaussian elimination. - ! A(N,N) is the original matrix in the input and transformed matrix - ! plus the pivoting element ratios below the diagonal in the output. - ! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001. - ! - IMPLICIT NONE - INTEGER, INTENT (IN) :: N - INTEGER :: I,J,K,ITMP - INTEGER, INTENT (OUT), DIMENSION (N) :: INDX - REAL (kind=RKIND) :: C1,PI,PI1,PJ - REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A - REAL (kind=RKIND), DIMENSION (N) :: C - ! - ! Initialize the index - ! - DO I = 1, N - INDX(I) = I - END DO - ! - ! Find the rescaling factors, one from each row - ! - DO I = 1, N - C1= 0.0 - DO J = 1, N - C1 = MAX(C1,ABS(A(I,J))) - END DO - C(I) = C1 - END DO - ! - ! Search the pivoting (largest) element from each column - ! - DO J = 1, N-1 - PI1 = 0.0 - DO I = J, N - PI = ABS(A(INDX(I),J))/C(INDX(I)) - IF (PI.GT.PI1) THEN - PI1 = PI - K = I - ENDIF - END DO - ! - ! Interchange the rows via INDX(N) to record pivoting order - ! - ITMP = INDX(J) - INDX(J) = INDX(K) - INDX(K) = ITMP - DO I = J+1, N - PJ = A(INDX(I),J)/A(INDX(J),J) - ! - ! Record pivoting ratios below the diagonal - ! - A(INDX(I),J) = PJ - ! - ! Modify other elements accordingly - ! - DO K = J+1, N - A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K) - END DO - END DO - END DO - ! - END SUBROUTINE ELGS - subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere_radius ) ! diff --git a/src/core_sw/mpas_sw_advection.F b/src/core_sw/mpas_sw_advection.F index b1ea06d385..5e04003375 100644 --- a/src/core_sw/mpas_sw_advection.F +++ b/src/core_sw/mpas_sw_advection.F @@ -12,6 +12,7 @@ module sw_advection use mpas_pool_routines use mpas_constants use mpas_log + use mpas_matrix_operations, only : mpas_migs contains @@ -635,14 +636,14 @@ subroutine sw_poly_fit_2(a_in,b_out,weights_in,m,n,ne) ata = matmul(at,a) ! if (m == n) then -! call sw_migs(a,n,b,indx) +! call mpas_migs(a,n,b,indx) ! else - call sw_migs(atha,n,atha_inv,indx) + call mpas_migs(atha,n,atha_inv,indx) b = matmul(atha_inv,ath) -! call sw_migs(ata,n,ata_inv,indx) +! call mpas_migs(ata,n,ata_inv,indx) ! b = matmul(ata_inv,at) ! end if b_out(1:n,1:m) = b(1:n,1:m) @@ -655,133 +656,6 @@ subroutine sw_poly_fit_2(a_in,b_out,weights_in,m,n,ne) end subroutine sw_poly_fit_2 - -! Updated 10/24/2001. -! -!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! ! -! Please Note: ! -! ! -! (1) This computer program is written by Tao Pang in conjunction with ! -! his book, "An Introduction to Computational Physics," published ! -! by Cambridge University Press in 1997. ! -! ! -! (2) No warranties, express or implied, are made for this program. ! -! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -subroutine sw_migs (A,N,X,INDX) -! -! subroutine to invert matrix A(N,N) with the inverse stored -! in X(N,N) in the output. Copyright (c) Tao Pang 2001. -! - IMPLICIT NONE - INTEGER, INTENT (IN) :: N - INTEGER :: I,J,K - INTEGER, INTENT (OUT), DIMENSION (N) :: INDX - REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A - REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X - REAL (kind=RKIND), DIMENSION (N,N) :: B -! - DO I = 1, N - DO J = 1, N - B(I,J) = 0.0 - END DO - END DO - DO I = 1, N - B(I,I) = 1.0 - END DO -! - call sw_elgs (A,N,INDX) -! - DO I = 1, N-1 - DO J = I+1, N - DO K = 1, N - B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K) - END DO - END DO - END DO -! - DO I = 1, N - X(N,I) = B(INDX(N),I)/A(INDX(N),N) - DO J = N-1, 1, -1 - X(J,I) = B(INDX(J),I) - DO K = J+1, N - X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I) - END DO - X(J,I) = X(J,I)/A(INDX(J),J) - END DO - END DO -end subroutine sw_migs - - -subroutine sw_elgs (A,N,INDX) -! -! subroutine to perform the partial-pivoting Gaussian elimination. -! A(N,N) is the original matrix in the input and transformed matrix -! plus the pivoting element ratios below the diagonal in the output. -! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001. -! - IMPLICIT NONE - INTEGER, INTENT (IN) :: N - INTEGER :: I,J,K,ITMP - INTEGER, INTENT (OUT), DIMENSION (N) :: INDX - REAL (kind=RKIND) :: C1,PI,PI1,PJ - REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A - REAL (kind=RKIND), DIMENSION (N) :: C -! -! Initialize the index -! - DO I = 1, N - INDX(I) = I - END DO -! -! Find the rescaling factors, one from each row -! - DO I = 1, N - C1= 0.0 - DO J = 1, N - C1 = MAX(C1,ABS(A(I,J))) - END DO - C(I) = C1 - END DO -! -! Search the pivoting (largest) element from each column -! - DO J = 1, N-1 - PI1 = 0.0 - DO I = J, N - PI = ABS(A(INDX(I),J))/C(INDX(I)) - IF (PI.GT.PI1) THEN - PI1 = PI - K = I - ENDIF - END DO -! -! Interchange the rows via INDX(N) to record pivoting order -! - ITMP = INDX(J) - INDX(J) = INDX(K) - INDX(K) = ITMP - DO I = J+1, N - PJ = A(INDX(I),J)/A(INDX(J),J) -! -! Record pivoting ratios below the diagonal -! - A(INDX(I),J) = PJ -! -! Modify other elements accordingly -! - DO K = J+1, N - A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K) - END DO - END DO - END DO -! -end subroutine sw_elgs - !------------------------------------------------------------- subroutine sw_initialize_deformation_weights( meshPool ) diff --git a/src/operators/Makefile b/src/operators/Makefile index 6a3d1301e7..2ad1c43384 100644 --- a/src/operators/Makefile +++ b/src/operators/Makefile @@ -21,7 +21,7 @@ operators: $(OBJS) $(DEPS) mpas_vector_operations.o: $(DEPS) mpas_matrix_operations.o: $(DEPS) mpas_tensor_operations.o: mpas_vector_operations.o mpas_matrix_operations.o $(DEPS) -mpas_rbf_interpolation.o: mpas_vector_operations.o +mpas_rbf_interpolation.o: mpas_vector_operations.o mpas_matrix_operations.o mpas_vector_reconstruction.o: mpas_rbf_interpolation.o mpas_spline_interpolation: mpas_tracer_advection_helpers.o: mpas_geometry_utils.o $(DEPS) diff --git a/src/operators/mpas_matrix_operations.F b/src/operators/mpas_matrix_operations.F index 953ee05f61..0f9e71c186 100644 --- a/src/operators/mpas_matrix_operations.F +++ b/src/operators/mpas_matrix_operations.F @@ -46,7 +46,7 @@ module mpas_matrix_operations mpas_matrix_cell_to_edge, & mpas_outer_product, & mpas_migs, & - mpas_elgs + mpas_legs !-------------------------------------------------------------------- ! @@ -437,132 +437,198 @@ subroutine mpas_outer_product(u,v,A)!{{{ end subroutine mpas_outer_product!}}} - ! Updated 10/24/2001. - ! - !!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ! - ! Please Note: ! - ! ! - ! (1) This computer program is written by Tao Pang in conjunction with ! - ! his book, "An Introduction to Computational Physics," published ! - ! by Cambridge University Press in 1997. ! - ! ! - ! (2) No warranties, express or implied, are made for this program. ! - ! ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - subroutine mpas_migs (a,n,x,indx)!{{{ - ! - ! subroutine to invert matrix a(n,n) with the inverse stored - ! in x(n,n) in the output. copyright (c) tao pang 2001. - ! - implicit none - integer, intent (in) :: n - integer :: i,j,k - integer, intent (out), dimension (n) :: indx - real (kind=RKIND), intent (inout), dimension (n,n):: a - real (kind=RKIND), intent (out), dimension (n,n):: x - real (kind=RKIND), dimension (n,n) :: b - ! - do i = 1, n - do j = 1, n - b(i,j) = 0.0 - end do - end do - do i = 1, n - b(i,i) = 1.0 - end do - ! - call mpas_elgs (a,n,indx) - ! - do i = 1, n-1 - do j = i+1, n - do k = 1, n - b(indx(j),k) = b(indx(j),k)-a(indx(j),i)*b(indx(i),k) - end do - end do - end do - ! - do i = 1, n - x(n,i) = b(indx(n),i)/a(indx(n),n) - do j = n-1, 1, -1 - x(j,i) = b(indx(j),i) - do k = j+1, n - x(j,i) = x(j,i)-a(indx(j),k)*x(k,i) - end do - x(j,i) = x(j,i)/a(indx(j),j) - end do - end do - end subroutine mpas_migs!}}} +!------------------------------------------------------------------------ +! +! The LU factorization and solve routines below (mpas_getrf, mpas_getrs) +! are adapted from the LAPACK reference implementation of DGETF2 and +! DGETRS, converted to free-form Fortran at MPAS real precision (RKIND). +! +! Copyright (c) 1992-2025 The University of Tennessee and The University +! of Tennessee Research Foundation. All rights +! reserved. +! Copyright (c) 2000-2025 The University of California Berkeley. All +! rights reserved. +! Copyright (c) 2006-2025 The University of Colorado Denver. All rights +! reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions are +! met: +! +! - Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! - Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions and the following disclaimer listed +! in this license in the documentation and/or other materials +! provided with the distribution. +! +! - Neither the name of the copyright holders nor the names of its +! contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! The copyright holders provide no reassurances that the source code +! provided does not infringe any patent, copyright, or any other +! intellectual property rights of third parties. The copyright holders +! disclaim any liability to any recipient for claims brought against +! recipient by any third party for infringement of that parties +! intellectual property rights. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +!------------------------------------------------------------------------ + + subroutine mpas_getrf(N, A, IPIV, INFO)!{{{ + + ! LU factorization with partial pivoting: A = P*L*U, in place. + + IMPLICIT NONE + + INTEGER, INTENT(IN):: N + REAL(KIND=RKIND), DIMENSION(N,N), INTENT(INOUT):: A + INTEGER, DIMENSION(N), INTENT(OUT):: IPIV + INTEGER, INTENT(OUT):: INFO + + REAL(KIND=RKIND), PARAMETER:: ZERO = 0.0_RKIND + + INTEGER:: I,J,JP + REAL(KIND=RKIND):: TMP + + INFO = 0 + + DO J = 1, N + ! Find pivot and test for singularity + + JP = J - 1 + MAXLOC(ABS(A(J:N,J)), DIM=1) + IPIV(J) = JP + IF(A(JP,J).NE.ZERO) THEN + ! Apply the interchange to columns 1:n + + IF (JP.NE.J) THEN + DO I = 1, N + TMP = A(J,I) + A(J,I) = A(JP,I) + A(JP,I) = TMP + ENDDO + ENDIF + + ! Compute elements j+1:n of j-th column + + IF (J.LT.N) THEN + DO I = 1, N-J + A(J+I,J) = A(J+I,J) / A(J,J) + ENDDO + ENDIF + + ELSE IF (INFO.EQ.0) THEN + INFO = J + ENDIF + + IF(J.LT.N) THEN + !Update trailing submatrix + + DO I = J+1, N + A(J+1:N,I) = A(J+1:N,I) - A(J+1:N,J) * A(J,I) + ENDDO + ENDIF + ENDDO + end subroutine mpas_getrf!}}} + + subroutine mpas_getrs(N, NRHS,A,IPIV,B)!{{{ + ! Solves A*X = B for nrhs right-hand sides using the P*L*U + ! factors computed by mpas_getrf + + IMPLICIT NONE + + INTEGER, INTENT(IN):: N,NRHS + REAL(KIND=RKIND), DIMENSION(N,N), INTENT(IN):: A + INTEGER, DIMENSION(N), INTENT(IN):: IPIV + REAL(KIND=RKIND), DIMENSION(N,NRHS), INTENT(INOUT):: B + + INTEGER:: I,K + REAL(KIND=RKIND):: TMP + + ! Apply row interchanges to the right hand sides + DO I = 1, N + IF (IPIV(I).NE.I) THEN + DO K = 1, NRHS + TMP = B(I,K) + B(I,K) = B(IPIV(I),K) + B(IPIV(I),K) = TMP + ENDDO + ENDIF + ENDDO + + ! Solve L*X = B, overwriting B with X + DO K = 1, NRHS + DO I =2, N + B(I,K) = B(I,K) - DOT_PRODUCT(A(I,1:I-1), B(1:I-1,K)) + ENDDO + ENDDO + + ! Solve U*X = B, overwriting B with X + DO K=1, NRHS + DO I = N, 1, -1 + B(I,K) = (B(I,K) - DOT_PRODUCT(A(I,I+1:N), B(I+1:N,K)))/A(I,I) + ENDDO + ENDDO + end subroutine mpas_getrs!}}} + + subroutine mpas_migs(A,N,X,INDX)!{{{ + ! Invert the N x N matrix A, returning the inverse in X + + IMPLICIT NONE + + INTEGER, INTENT(IN):: N + REAL(KIND=RKIND), DIMENSION(N,N), INTENT(INOUT):: A + REAL(KIND=RKIND), DIMENSION(N,N), INTENT(OUT):: X + INTEGER, DIMENSION(N), INTENT(OUT):: INDX + + INTEGER::K, INFO + + CALL mpas_getrf(N,A,INDX,INFO) + + X(:,:) = 0.0_RKIND + DO K = 1,N + X(K,K) = 1.0_RKIND + ENDDO + + CALL mpas_getrs(N,N,A,INDX,X) - subroutine mpas_elgs (a,n,indx)!{{{ - ! - ! subroutine to perform the partial-pivoting gaussian elimination. - ! a(n,n) is the original matrix in the input and transformed matrix - ! plus the pivoting element ratios below the diagonal in the output. - ! indx(n) records the pivoting order. copyright (c) tao pang 2001. - ! - implicit none - integer, intent (in) :: n - integer :: i,j,k,itmp - integer, intent (out), dimension (n) :: indx - real (kind=RKIND) :: c1,pi,pi1,pj - real (kind=RKIND), intent (inout), dimension (n,n) :: a - real (kind=RKIND), dimension (n) :: c - ! - ! initialize the index - ! - do i = 1, n - indx(i) = i - end do - ! - ! find the rescaling factors, one from each row - ! - do i = 1, n - c1= 0.0 - do j = 1, n - c1 = max(c1,abs(a(i,j))) - end do - c(i) = c1 - end do - ! - ! search the pivoting (largest) element from each column - ! - do j = 1, n-1 - pi1 = 0.0 - do i = j, n - pi = abs(a(indx(i),j))/c(indx(i)) - if (pi.gt.pi1) then - pi1 = pi - k = i - endif - end do - ! - ! interchange the rows via indx(n) to record pivoting order - ! - itmp = indx(j) - indx(j) = indx(k) - indx(k) = itmp - do i = j+1, n - pj = a(indx(i),j)/a(indx(j),j) - ! - ! record pivoting ratios below the diagonal - ! - a(indx(i),j) = pj - ! - ! modify other elements accordingly - ! - do k = j+1, n - a(indx(i),k) = a(indx(i),k)-pj*a(indx(j),k) - end do - end do - end do - ! - end subroutine mpas_elgs!}}} + end subroutine mpas_migs!}}} + subroutine mpas_legs(A,N,B,X,INDX)!{{{ + ! Solve the linear system A*X = B by LU factorization. + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: N + REAL(KIND=RKIND), DIMENSION(N,N), INTENT(INOUT) :: A + REAL(KIND=RKIND), DIMENSION(N), INTENT(INOUT) :: B + REAL(KIND=RKIND), DIMENSION(N), INTENT(OUT) :: X + INTEGER, DIMENSION(N), INTENT(OUT) :: INDX + + REAL(KIND=RKIND), DIMENSION(N,1) :: RHS + INTEGER :: INFO + + CALL mpas_getrf(N,A,INDX,INFO) + + RHS(:,1) = B(:) + CALL mpas_getrs(N,1,A,INDX,RHS) + X(:) = RHS(:,1) + + end subroutine mpas_legs!}}} end module mpas_matrix_operations diff --git a/src/operators/mpas_rbf_interpolation.F b/src/operators/mpas_rbf_interpolation.F index 1c136a3c91..3b42919b7f 100644 --- a/src/operators/mpas_rbf_interpolation.F +++ b/src/operators/mpas_rbf_interpolation.F @@ -21,6 +21,7 @@ module mpas_rbf_interpolation use mpas_dmpar use mpas_derived_types use mpas_vector_operations + use mpas_matrix_operations, only : mpas_legs implicit none private @@ -1616,234 +1617,5 @@ subroutine mpas_set_up_vector_free_slip_rbf_matrix_and_rhs(pointCount, dimension end subroutine mpas_set_up_vector_free_slip_rbf_matrix_and_rhs!}}} -! Updated 10/24/2001. -! -!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! ! -! Please Note: ! -! ! -! (1) This computer program is written by Tao Pang in conjunction with ! -! his book, "An Introduction to Computational Physics," published ! -! by Cambridge University Press in 1997. ! -! ! -! (2) No warranties, express or implied, are made for this program. ! -! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!PROGRAM EX43 -! -! -! An example of solving linear equation set A(N,N)*X(N) = B(N) -! with the partial-pivoting Gaussian elimination scheme. The -! numerical values are for the Wheatstone bridge example discussed -! in Section 4.1 in the book with all resistances being 100 ohms -! and the voltage 200 volts. -! -! IMPLICIT NONE -! INTEGER, PARAMETER :: N=3 -! INTEGER :: I,J -! INTEGER, DIMENSION (N) :: INDX -! REAL, DIMENSION (N) :: X,B -! REAL, DIMENSION (N,N) :: A -! DATA B /200.0,0.0,0.0/, & -! ((A(I,J), J=1,N),I=1,N) /100.0,100.0,100.0,-100.0, & -! 300.0,-100.0,-100.0,-100.0, 300.0/ -! -! call mpas_legs (A,N,B,X,INDX) -! -! WRITE (6, "(F16.8)") (X(I), I=1,N) -!END PROGRAM EX43 - -!*********************************************************************** -! -! routine mpas_legs -! -!> \brief MPAS Gaussian elimination solver routine -!> \author Xylar Asay-Davis -!> \date 03/28/13 -!> \details -!> This routine solves the equation A(N,N)*X(N) = B(N) with the partial-pivoting -!> Gaussian Elimination scheme. Copyright (c) Tao Pang 2001. -!----------------------------------------------------------------------- -subroutine mpas_legs (A,N,B,X,INDX)!{{{ - - IMPLICIT NONE - integer, INTENT (IN) :: N !< Input: Size of matrix and vectors - integer, INTENT (OUT), DIMENSION (N) :: INDX !< Output: Pivot vector - real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A !< Input/Output: Matrix - real(kind=RKIND), INTENT (INOUT), DIMENSION (N) :: B !< Input/Output: Right hand side vector - real(kind=RKIND), INTENT (OUT), DIMENSION (N) :: X !< Output: Solution - - integer :: I,J -! - CALL elgs (A,N,INDX) -! - DO I = 1, N-1 - DO J = I+1, N - B(INDX(J)) = B(INDX(J))-A(INDX(J),I)*B(INDX(I)) - END DO - END DO -! - X(N) = B(INDX(N))/A(INDX(N),N) - DO I = N-1, 1, -1 - X(I) = B(INDX(I)) - DO J = I+1, N - X(I) = X(I)-A(INDX(I),J)*X(J) - END DO - X(I) = X(I)/A(INDX(I),I) - END DO -! -END subroutine mpas_legs!}}} -! - - - -! Updated 10/24/2001. -! -!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! ! -! Please Note: ! -! ! -! (1) This computer program is written by Tao Pang in conjunction with ! -! his book, "An Introduction to Computational Physics," published ! -! by Cambridge University Press in 1997. ! -! ! -! (2) No warranties, express or implied, are made for this program. ! -! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!*********************************************************************** -! -! routine migs -! -!> \brief Matrix inversion routine -!> \author Xylar Asay-Davis -!> \date 03/28/13 -!> \details -!> This routine inverts the matrix A(N,N) and stores it in X(N,B) -!> Copyright (c) Tao Pang 2001. -!----------------------------------------------------------------------- -subroutine migs (A,N,X,INDX)!{{{ - IMPLICIT NONE - integer, INTENT (IN) :: N !< Input: Size of matrix and inverse - integer, INTENT (OUT), DIMENSION (N) :: INDX !< Output: Pivot vector - real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A !< Input/Output: Matrix to invert - real(kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X !< Output: Inverse of Matrix - real(kind=RKIND), DIMENSION (N,N) :: B - integer :: I,J,K -! - DO I = 1, N - DO J = 1, N - B(I,J) = 0.0 - END DO - END DO - DO I = 1, N - B(I,I) = 1.0 - END DO -! - CALL elgs (A,N,INDX) -! - DO I = 1, N-1 - DO J = I+1, N - DO K = 1, N - B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K) - END DO - END DO - END DO -! - DO I = 1, N - X(N,I) = B(INDX(N),I)/A(INDX(N),N) - DO J = N-1, 1, -1 - X(J,I) = B(INDX(J),I) - DO K = J+1, N - X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I) - END DO - X(J,I) = X(J,I)/A(INDX(J),J) - END DO - END DO -END subroutine migs!}}} - -!*********************************************************************** -! -! routine elgs -! -!> \brief Partial-pivoting Gaussian elimination routine -!> \author Xylar Asay-Davis -!> \date 03/28/13 -!> \details -!> This routine performs the partial-pivoting Gaussian elimination. -!> Copyright (c) Tao Pang 2001. -!----------------------------------------------------------------------- - -subroutine elgs (A,N,INDX)!{{{ -! -! subroutine to perform the partial-pivoting Gaussian elimination. -! A(N,N) is the original matrix in the input and transformed matrix -! plus the pivoting element ratios below the diagonal in the output. -! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001. -! - IMPLICIT NONE - integer, INTENT (IN) :: N !< Input: Size of matrix - integer, INTENT (OUT), DIMENSION (N) :: INDX !< Output: Pivot vector - real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A !< Input/Output: Matrix and solution - integer :: I,J,K,ITMP - real(kind=RKIND) :: C1,PI,PI1,PJ - real(kind=RKIND), DIMENSION (N) :: C -! -! Initialize the index -! - DO I = 1, N - INDX(I) = I - END DO -! -! Find the rescaling factors, one from each row -! - DO I = 1, N - C1= 0.0 - DO J = 1, N - !C1 = AMAX1(C1,ABS(A(I,J))) - C1 = MAX(C1,ABS(A(I,J))) - END DO - C(I) = C1 - END DO -! -! Search the pivoting (largest) element from each column -! - DO J = 1, N-1 - PI1 = 0.0 - DO I = J, N - PI = ABS(A(INDX(I),J))/C(INDX(I)) - IF (PI.GT.PI1) THEN - PI1 = PI - K = I - ENDIF - END DO -! -! Interchange the rows via INDX(N) to record pivoting order -! - ITMP = INDX(J) - INDX(J) = INDX(K) - INDX(K) = ITMP - DO I = J+1, N - PJ = A(INDX(I),J)/A(INDX(J),J) -! -! Record pivoting ratios below the diagonal -! - A(INDX(I),J) = PJ -! -! Modify other elements accordingly -! - DO K = J+1, N - A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K) - END DO - END DO - END DO -! -END subroutine elgs!}}} - end module mpas_rbf_interpolation