Skip to content

Replaced linear system matrix operation functions with modified LAPACK implementations. #1487

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
robertmullis wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Replaced linear system matrix operation functions with modified LAPACK implementations. #1487

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a3da0e7
Refactored matrix operations using LAPACK library
Jul 9, 2026
a12a127
added copyright disclaimer to the mpas_getrf and mpas_getrs functions
Jul 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
131 changes: 2 additions & 129 deletions src/core_init_atmosphere/mpas_atm_advection.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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 )

!
Expand Down
134 changes: 4 additions & 130 deletions src/core_sw/mpas_sw_advection.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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 )
Expand Down
2 changes: 1 addition & 1 deletion src/operators/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading