Merge pull request #58 from hphratchian/master

Added functionality for reordering matrix columns in MQC_General and …

src/mqc_algebra2.F03

@@ -59,7 +59,6 @@ Module MQC_Algebra2
         integer(kind=int64),private::rank=-1
         character(len=64),private::dataType,storageFormat
         integer(kind=int64),dimension(10),private::dimensions=0
-!hph        real(kind=real64),dimension(:),allocatable,private::realArray
         real(kind=real64),dimension(:),allocatable::realArray
         integer(kind=int64),dimension(:),allocatable,private::integerArray
         logical,private::initialized=.false.
@@ -1128,6 +1127,65 @@ subroutine MQC_Variable_MatrixPermuteColumns(mqcVariable,columnI,columnJ)
       end subroutine MQC_Variable_MatrixPermuteColumns
 
 
+!hph+
+!!
+!!PROCEDURE MQC_Variable_MatrixOrderedColumns
+!      subroutine MQC_Variable_MatrixOrderedColumns(mqcVariable,map)
+!!
+!!     This subroutine is used to order the columns of a matrix according to the
+!!     array map. If this routine is called with an mqcVariable that isn't rank
+!!     2, it results in error.
+!!
+!!
+!!     H. P. Hratchian, 2023.
+!!
+!!
+!!     Variable Declarations.
+!      implicit none
+!      class(MQC_Variable)::mqcVariable
+!      integer(kind=int64),intent(in)::columnI,columnJ
+!      integer(kind=int64)::lenVector,iStart,jStart,iEnd,jEnd
+!      integer(kind=int64),dimension(:),allocatable::tmpArrayI
+!      real(kind=real64),dimension(:),allocatable::tmpArrayR
+!      integer(kind=int64)::k
+!!
+!!
+!!     Do the work...
+!!
+!      if(.not.mqcVariable%initialized) call mqc_error('Cannot use MQC_Variable_Get with uninitialized permuteColumns.')
+!      if(mqcVariable%rank.ne.2)  &
+!        call mqc_error('permuteColumns: Rank of variable must be 2.')
+!      if(TRIM(mqcVariable%storageFormat).ne.'FULL')   &
+!        call mqc_error_a('permuteColumns: Can only be used on FULL storage arrays.',  &
+!        6,'Current storage form:',TRIM(mqcVariable%storageFormat))
+!      lenVector = SIZE(mqcVariable,1)
+!      iStart = MQC_Variable_getArrayPosition(mqcVariable,[ 1, columnI ])
+!      jStart = MQC_Variable_getArrayPosition(mqcVariable,[ 1, columnJ ])
+!      iEnd = iStart + lenVector - 1
+!      jEnd = jStart + lenVector - 1
+!      select case(TRIM(MQC_Variable_getType(mqcVariable)))
+!      case('REAL')
+!        Allocate(tmpArrayR(lenVector))
+!        tmpArrayR = mqcVariable%realArray(iStart:iEnd)
+!        mqcVariable%realArray(iStart:iEnd) = mqcVariable%realArray(jStart:jEnd)
+!        mqcVariable%realArray(jStart:jEnd) = tmpArrayR
+!        DeAllocate(tmpArrayR)
+!      case('INTEGER')
+!        Allocate(tmpArrayI(lenVector))
+!        tmpArrayI = mqcVariable%integerArray(iStart:iEnd)
+!        mqcVariable%integerArray(iStart:iEnd) = mqcVariable%integerArray(jStart:jEnd)
+!        mqcVariable%integerArray(jStart:jEnd) = tmpArrayI
+!        DeAllocate(tmpArrayI)
+!      case default
+!        call mqc_error_a('permuteColumns: Unknown MQC variable type.', 6, &
+!          'TRIM(MQC_Variable_getType(mqcVariable))', TRIM(MQC_Variable_getType(mqcVariable)) )
+!      end select
+!!
+!      return
+!      end subroutine MQC_Variable_MatrixPermuteColumns
+!hph-
+
+
 !
 !PROCEDURE MQC_Variable_put_MQC
       subroutine MQC_Variable_put_MQC(mqcVariable,valueMQC,arrayElement)
@@ -1664,10 +1722,6 @@ subroutine MQC_Variable_setVal(mqcVariable,scalarIntegerIn,  &
             myDimensions(1) = SIZE(arrayIntegerIn)
           elseIf(fillReal) then
             myDimensions(1) = SIZE(arrayRealIn)
-
-            write(*,*)'    myRank = ',myRank
-            write(*,*)'    myDimensions(1) = ',myDimensions(1)
-
           else
             call mqc_error_l(  &
               'MQC_Variable_setVal: Confused setting myDimensions.', 6, &
@@ -1683,7 +1737,6 @@ subroutine MQC_Variable_setVal(mqcVariable,scalarIntegerIn,  &
             'fillArray', fillArray )
         endIf
       endIf
-      if(DEBUG) write(*,*)'   myRank = ',myRank
 !
 !     Initialize mqcVariable and then set the value appropriately.
 !
@@ -2081,7 +2134,7 @@ subroutine MQC_Variable_mqc2intrinsicInteger1Array(intrinsicOut,mqcVariable)
       case(3)
         intrinsicOut = mqcVariable%integerArray(:)
       case default
-        call mqc_error_I('MQC_Variable_mqc2intrinsicReal1Array: '// &
+        call mqc_error_I('MQC_Variable_mqc2intrinsicInteger1Array: '// &
           'Unknown MQC_Variable type found.',6, &
           'MQC_Variable_getTypeCode(mqcVariable)',  &
           MQC_Variable_getTypeCode(mqcVariable))

src/mqc_gaussian.F03

@@ -791,7 +791,7 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Header(fileinfo,  &
           if(.not.ok) Call MQC_Error_L('Error opening Gaussian matrix file.', 6, &
              'ok', ok )
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (1): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -919,7 +919,7 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Write_Header(fileinfo,  &
           if(.not.ok) Call MQC_Error_L('Error opening Gaussian matrix file.', 6, &
                'ok', ok )
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (2): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -1152,7 +1152,7 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Read_Array(fileinfo,  &
           call MQC_Gaussian_Unformatted_Matrix_Read_Header(fileinfo,  &
             filename)
         else
-          call MQC_Error_l('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_l('Error reading Gaussian matrix file header (3): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -1648,7 +1648,7 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Write_Array2(fileinfo,  &
                'ok', ok )
           call MQC_Gaussian_Unformatted_Matrix_Write_Header(fileinfo,filename)
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (4): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -2080,7 +2080,7 @@ subroutine MQC_Gaussian_Unformatted_Matrix_Write_Array(fileinfo,  &
             'ok', ok )
           call MQC_Gaussian_Unformatted_Matrix_Write_Header(fileinfo,filename)
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (5): Must include a filename.', 6, &
             'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -3486,7 +3486,7 @@ subroutine mqc_gaussian_unformatted_matrix_get_EST_object(fileinfo,label, &
           call MQC_Gaussian_Unformatted_Matrix_Read_Header(fileinfo,  &
             filename)
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (6): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -5145,7 +5145,7 @@ subroutine mqc_gaussian_unformatted_matrix_get_twoERIs(fileinfo,label, &
           call MQC_Gaussian_Unformatted_Matrix_Read_Header(fileinfo,  &
             filename)
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (7): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf
@@ -5363,7 +5363,7 @@ subroutine mqc_gaussian_unformatted_matrix_write_twoERIs(fileinfo,label, &
                'ok', ok )
           call MQC_Gaussian_Unformatted_Matrix_Write_Header(fileinfo,filename)
         else
-          call MQC_Error_L('Error reading Gaussian matrix file header: Must include a filename.', 6, &
+          call MQC_Error_L('Error reading Gaussian matrix file header (8): Must include a filename.', 6, &
                'PRESENT(filename)', PRESENT(filename) )
         endIf
       endIf

src/mqc_general.F03

@@ -133,6 +133,11 @@ Module MQC_General
         module procedure mqc_bubbleSort_integer
         module procedure mqc_bubbleSort_real
       end interface
+!
+      interface matrixOrderedColumns
+        module procedure mqc_matrixOrderedColumns_integer
+        module procedure mqc_matrixOrderedColumns_real
+      end interface
 !
       interface flatten
         module procedure mqc_flattenR4Real
@@ -2055,7 +2060,7 @@ end subroutine mqc_matrixInverse_symmFull
 !----------------------------------------------------------------
 !
 !PROCEDURE mqc_bubbleSort_integer
-      subroutine mqc_bubbleSort_integer(listIn,listOut,map,sortListIn)
+      subroutine mqc_bubbleSort_integer(listIn,listOut,map,sortListIn,reverse)
 !
 !     This subroutine carries out a simple bubble sort algorithm to order the
 !     values in the array <listIn>. The sorted list is returned in the optional
@@ -2064,7 +2069,9 @@ subroutine mqc_bubbleSort_integer(listIn,listOut,map,sortListIn)
 !     <listOut> or <map> is sent, then the sorted list will overwrite the input
 !     data in listIn. If <listOut> and/or <map> are sent, it should already be
 !     appropriately allocated. If <sortListIn> is sent and is TRUE, then listIn
-!     is always returned with the sorted values.
+!     is always returned with the sorted values. Optional argument <reverse> is
+!     sent as TRUE if the list should be sorted from largest to smallest value.
+!     The default value of <reverse> is FALSE.
 !
 !     H. P. Hratchian, 2023.
 !
@@ -2073,18 +2080,25 @@ subroutine mqc_bubbleSort_integer(listIn,listOut,map,sortListIn)
       implicit none
       integer(kind=int64),dimension(:)::listIn
       integer(kind=int64),dimension(:),optional::listOut,map
-      logical,optional::sortListIn
+      logical,optional::sortListIn,reverse
       integer(kind=int64)::i,j,nDim,nSwaps,valueTemp
       integer(kind=int64),dimension(:),allocatable::listTemp,listMapTemp
-      logical::overwriteListIn
+      logical::overwriteListIn,doReverse
 !
-!     Allocate listTemp and copy listIn into it.
+!     Set up overwriteListIn and doReverse based on whether the related optional
+!     dummy arguments have been sent by the calling program unit.
 !
       overwriteListIn = .false.
       if(PRESENT(sortListIn)) overwriteListIn = sortListIn
+      doReverse = .false.
+      if(PRESENT(reverse)) doReverse = reverse
+!
+!     Allocate listTemp and copy listIn into it.
+!
       nDim = Size(listIn)
       Allocate(listTemp(nDim),listMapTemp(nDim))
       listTemp = listIn
+      if(doReverse) listTemp = -listTemp
       call mqc_seq(listMapTemp)
 !
 !     Carry out the bubble sort algorithm on listTemp.
@@ -2107,6 +2121,7 @@ subroutine mqc_bubbleSort_integer(listIn,listOut,map,sortListIn)
 !
 !     Put the sorted list into listOut or back into listIn.
 !
+      if(doReverse) listTemp = -listTemp
       if(PRESENT(listOut)) listOut = listTemp
       if(PRESENT(map)) map = listMapTemp
       if(overwriteListIn.or.(.not.(PRESENT(listOut).or.PRESENT(map))))  &
@@ -2118,7 +2133,7 @@ end subroutine mqc_bubbleSort_integer
 
 !
 !PROCEDURE mqc_bubbleSort_real
-      subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn)
+      subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn,reverse)
 !
 !     This subroutine carries out a simple bubble sort algorithm to order the
 !     values in the array <listIn>. The sorted list is returned in the optional
@@ -2127,7 +2142,10 @@ subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn)
 !     <listOut> or <map> is sent, then the sorted list will overwrite the input
 !     data in listIn. If <listOut> and/or <map> are sent, it should already be
 !     appropriately allocated. If <sortListIn> is sent and is TRUE, then listIn
-!     is always returned with the sorted values.
+!     is always returned with the sorted values. Optional argument <reverse> is
+!     sent as TRUE if the list should be sorted from largest to smallest value.
+!     The default value of <reverse> is FALSE.
+!
 !
 !     H. P. Hratchian, 2023.
 !
@@ -2137,20 +2155,27 @@ subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn)
       real(kind=real64),dimension(:)::listIn
       real(kind=real64),dimension(:),optional::listOut
       integer(kind=int64),dimension(:),optional::map
-      logical,optional::sortListIn
+      logical,optional::sortListIn,reverse
       integer(kind=int64)::i,j,nDim,nSwaps
       integer(kind=int64),dimension(:),allocatable::listMapTemp
       real(kind=real64)::valueTemp
       real(kind=real64),dimension(:),allocatable::listTemp
-      logical::overwriteListIn
+      logical::overwriteListIn,doReverse
 !
-!     Allocate listTemp and copy listIn into it.
+!     Set up overwriteListIn and doReverse based on whether the related optional
+!     dummy arguments have been sent by the calling program unit.
 !
       overwriteListIn = .false.
       if(PRESENT(sortListIn)) overwriteListIn = sortListIn
+      doReverse = .false.
+      if(PRESENT(reverse)) doReverse = reverse
+!
+!     Allocate listTemp and copy listIn into it.
+!
       nDim = Size(listIn)
       Allocate(listTemp(nDim),listMapTemp(nDim))
       listTemp = listIn
+      if(doReverse) listTemp = -listTemp
       call mqc_seq(listMapTemp)
 !
 !     Carry out the bubble sort algorithm on listTemp.
@@ -2173,6 +2198,7 @@ subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn)
 !
 !     Put the sorted list into listOut or back into listIn.
 !
+      if(doReverse) listTemp = -listTemp
       if(PRESENT(listOut)) listOut = listTemp
       if(PRESENT(map)) map = listMapTemp
       if(overwriteListIn.or.(.not.(PRESENT(listOut).or.PRESENT(map))))  &
@@ -2182,6 +2208,80 @@ subroutine mqc_bubbleSort_real(listIn,listOut,map,sortListIn)
       end subroutine mqc_bubbleSort_real
 
 
+!
+!PROCEDURE mqc_matrixOrderedColumns_integer
+      subroutine mqc_matrixOrderedColumns_integer(matrix,map)
+!
+!     This routine reordered the columns of the matrix <matrix> according to the
+!     mapping array <map>.
+!
+!
+!     H. P. Hratchian, 2023.
+!
+!
+!     Variable Declarations.
+      implicit none
+      integer(kind=int64),dimension(:,:)::matrix
+      integer(kind=int64),dimension(:)::map
+      integer(kind=int64)::i
+      integer(kind=int64),dimension(:,:),allocatable::tempMatrix
+
+!
+!     Begin by ensuring map is the same length as the number of columns in
+!     matrix.
+!
+      if(SIZE(matrix,2).ne.SIZE(map))  &
+        call mqc_error('mqc_matrixOrderedColumns: map is wrong length!')
+!
+!     Do the work.
+!
+      Allocate(tempMatrix(Size(matrix,1),Size(matrix,2)))
+      do i = 1,Size(matrix,2)
+        tempMatrix(:,i) = matrix(:,map(i))
+      endDo
+      matrix = tempMatrix
+!
+      return
+      end subroutine mqc_matrixOrderedColumns_integer
+
+
+!
+!PROCEDURE mqc_matrixOrderedColumns_real
+      subroutine mqc_matrixOrderedColumns_real(matrix,map)
+!
+!     This routine reordered the columns of the matrix <matrix> according to the
+!     mapping array <map>.
+!
+!
+!     H. P. Hratchian, 2023.
+!
+!
+!     Variable Declarations.
+      implicit none
+      real(kind=real64),dimension(:,:)::matrix
+      integer(kind=int64),dimension(:)::map
+      integer(kind=int64)::i
+      real(kind=real64),dimension(:,:),allocatable::tempMatrix
+
+!
+!     Begin by ensuring map is the same length as the number of columns in
+!     matrix.
+!
+      if(SIZE(matrix,2).ne.SIZE(map))  &
+        call mqc_error('mqc_matrixOrderedColumns: map is wrong length!')
+!
+!     Do the work.
+!
+      Allocate(tempMatrix(Size(matrix,1),Size(matrix,2)))
+      do i = 1,Size(matrix,2)
+        tempMatrix(:,i) = matrix(:,map(i))
+      endDo
+      matrix = tempMatrix
+!
+      return
+      end subroutine mqc_matrixOrderedColumns_real
+
+
 !
 !PROCEDURE mqc_flatten
       function mqc_flattenR4Real(inArray) result(outArray)