matrix_vector.f90

Go to the documentation of this file.

MODULE matrix_vector

  USE base_routines
  USE constants
  USE input_output
  USE iso_varying_string
  USE kinds
  USE lists
  USE strings
  USE types
  USE linkedlist_routines

#include "macros.h"  

  IMPLICIT NONE

  !PRIVATE

  !Module parameters

  INTEGER(INTG), PARAMETER :: matrix_vector_intg_type=integer_type 
  INTEGER(INTG), PARAMETER :: matrix_vector_sp_type=single_real_type 
  INTEGER(INTG), PARAMETER :: matrix_vector_dp_type=double_real_type 
  INTEGER(INTG), PARAMETER :: matrix_vector_l_type=logical_type 
  
  INTEGER(INTG), PARAMETER :: matrix_block_storage_type=0
  INTEGER(INTG), PARAMETER :: matrix_diagonal_storage_type=1
  INTEGER(INTG), PARAMETER :: matrix_column_major_storage_type=2
  INTEGER(INTG), PARAMETER :: matrix_row_major_storage_type=3
  INTEGER(INTG), PARAMETER :: matrix_compressed_row_storage_type=4
  INTEGER(INTG), PARAMETER :: matrix_compressed_column_storage_type=5
  INTEGER(INTG), PARAMETER :: matrix_row_column_storage_type=6

  INTEGER(INTG), PARAMETER :: bisectiontolinearsearchthreshold=10

  !Module types

  !Matrix types
  
  !Module variables

  INTEGER(INTG), SAVE :: matrix_vector_id=1

  !Interfaces

  INTERFACE matrix_all_values_set
    MODULE PROCEDURE matrix_all_values_set_intg
    MODULE PROCEDURE matrix_all_values_set_sp
    MODULE PROCEDURE matrix_all_values_set_dp
    MODULE PROCEDURE matrix_all_values_set_l
  END INTERFACE !MATRIX_ALL_VALUES_SET

  INTERFACE matrix_data_get
    MODULE PROCEDURE matrix_data_get_intg
    MODULE PROCEDURE matrix_data_get_sp
    MODULE PROCEDURE matrix_data_get_dp
    MODULE PROCEDURE matrix_data_get_l
  END INTERFACE !MATRIX_DATA_GET

  INTERFACE matrix_values_add
    MODULE PROCEDURE matrix_values_add_intg
    MODULE PROCEDURE matrix_values_add_intg1
    MODULE PROCEDURE matrix_values_add_intg2
    MODULE PROCEDURE matrix_values_add_sp
    MODULE PROCEDURE matrix_values_add_sp1
    MODULE PROCEDURE matrix_values_add_sp2
    MODULE PROCEDURE matrix_values_add_dp
    MODULE PROCEDURE matrix_values_add_dp1
    MODULE PROCEDURE matrix_values_add_dp2
    MODULE PROCEDURE matrix_values_add_l
    MODULE PROCEDURE matrix_values_add_l1
    MODULE PROCEDURE matrix_values_add_l2
  END INTERFACE !MATRIX_VALUES_ADD

  INTERFACE matrix_values_get
    MODULE PROCEDURE matrix_values_get_intg
    MODULE PROCEDURE matrix_values_get_intg1
    MODULE PROCEDURE matrix_values_get_intg2
    MODULE PROCEDURE matrix_values_get_sp
    MODULE PROCEDURE matrix_values_get_sp1
    MODULE PROCEDURE matrix_values_get_sp2
    MODULE PROCEDURE matrix_values_get_dp
    MODULE PROCEDURE matrix_values_get_dp1
    MODULE PROCEDURE matrix_values_get_dp2
    MODULE PROCEDURE matrix_values_get_l
    MODULE PROCEDURE matrix_values_get_l1
    MODULE PROCEDURE matrix_values_get_l2
  END INTERFACE !MATRIX_VALUES_GET
  
  INTERFACE matrix_values_set
    MODULE PROCEDURE matrix_values_set_intg
    MODULE PROCEDURE matrix_values_set_intg1
    MODULE PROCEDURE matrix_values_set_intg2
    MODULE PROCEDURE matrix_values_set_sp
    MODULE PROCEDURE matrix_values_set_sp1
    MODULE PROCEDURE matrix_values_set_sp2
    MODULE PROCEDURE matrix_values_set_dp
    MODULE PROCEDURE matrix_values_set_dp1
    MODULE PROCEDURE matrix_values_set_dp2
    MODULE PROCEDURE matrix_values_set_l
    MODULE PROCEDURE matrix_values_set_l1
    MODULE PROCEDURE matrix_values_set_l2
  END INTERFACE !MATRIX_VALUES_SET

  INTERFACE vector_all_values_set
    MODULE PROCEDURE vector_all_values_set_intg
    MODULE PROCEDURE vector_all_values_set_sp
    MODULE PROCEDURE vector_all_values_set_dp
    MODULE PROCEDURE vector_all_values_set_l
  END INTERFACE !VECTOR_ALL_VALUES_SET

  INTERFACE vector_data_get
    MODULE PROCEDURE vector_data_get_intg
    MODULE PROCEDURE vector_data_get_sp
    MODULE PROCEDURE vector_data_get_dp
    MODULE PROCEDURE vector_data_get_l
  END INTERFACE !VECTOR_DATA_GET

   INTERFACE vector_values_get
    MODULE PROCEDURE vector_values_get_intg
    MODULE PROCEDURE vector_values_get_intg1
    MODULE PROCEDURE vector_values_get_sp
    MODULE PROCEDURE vector_values_get_sp1
    MODULE PROCEDURE vector_values_get_dp
    MODULE PROCEDURE vector_values_get_dp1
    MODULE PROCEDURE vector_values_get_l
    MODULE PROCEDURE vector_values_get_l1
  END INTERFACE !VECTOR_VALUES_GET
  
  INTERFACE vector_values_set
    MODULE PROCEDURE vector_values_set_intg
    MODULE PROCEDURE vector_values_set_intg1
    MODULE PROCEDURE vector_values_set_sp
    MODULE PROCEDURE vector_values_set_sp1
    MODULE PROCEDURE vector_values_set_dp
    MODULE PROCEDURE vector_values_set_dp1
    MODULE PROCEDURE vector_values_set_l
    MODULE PROCEDURE vector_values_set_l1
  END INTERFACE !VECTOR_VALUES_SET

  PUBLIC matrix_vector_intg_type,matrix_vector_sp_type,matrix_vector_dp_type,matrix_vector_l_type

  PUBLIC matrix_block_storage_type,matrix_diagonal_storage_type,matrix_column_major_storage_type,matrix_row_major_storage_type, &
    & matrix_compressed_row_storage_type,matrix_compressed_column_storage_type,matrix_row_column_storage_type

  PUBLIC matrix_create_finish,matrix_create_start,matrix_data_get,matrix_datatypeget,matrix_data_type_set,matrix_destroy, &
    & matrix_duplicate,matrix_max_columns_per_row_get,matrix_number_non_zeros_set,matrix_number_non_zeros_get,matrix_max_size_set, &
    & matrix_output,matrix_size_set,matrix_storage_location_find,matrix_storage_locations_set,matrix_storage_type_get, &
    & matrix_storage_type_set,matrix_values_add,matrix_values_get,matrix_values_set

  PUBLIC vector_all_values_set,vector_create_finish,vector_create_start,vector_data_get,vector_datatypeget,vector_data_type_set, &
    & vector_destroy,vector_duplicate,vector_size_set,vector_values_get,vector_values_set

  PUBLIC matrix_linklist_set,matrix_linklist_get
CONTAINS
  
  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_all_values_set_intg(MATRIX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_ALL_VALUES_SET_INTG",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
          matrix%DATA_INTG=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_ALL_VALUES_SET_INTG")
    RETURN
999 errorsexits("MATRIX_ALL_VALUES_SET_INTG",err,error)
    RETURN 1
  END SUBROUTINE matrix_all_values_set_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_all_values_set_sp(MATRIX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    REAL(SP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_ALL_VALUES_SET_SP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
          matrix%DATA_SP=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_ALL_VALUES_SET_SP")
    RETURN
999 errorsexits("MATRIX_ALL_VALUES_SET_SP",err,error)
    RETURN 1
  END SUBROUTINE matrix_all_values_set_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_all_values_set_dp(MATRIX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_ALL_VALUES_SET_DP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
          matrix%DATA_DP=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_ALL_VALUES_SET_DP")
    RETURN
999 errorsexits("MATRIX_ALL_VALUES_SET_DP",err,error)
    RETURN 1
  END SUBROUTINE matrix_all_values_set_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_all_values_set_l(MATRIX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    LOGICAL, INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_ALL_VALUES_SET_L",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
          matrix%DATA_L=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_ALL_VALUES_SET_L")
    RETURN
999 errorsexits("MATRIX_ALL_VALUES_SET_L",err,error)
    RETURN 1
  END SUBROUTINE matrix_all_values_set_l

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_create_finish(MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: column_idx,COUNT,row_idx,row_idx2
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_CREATE_FINISH",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("Matrix has been finished.",err,error,*999)
      ELSE
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          IF(matrix%MAX_M==-1) matrix%MAX_M=matrix%M
          IF(matrix%MAX_N==-1) matrix%MAX_N=matrix%N
          matrix%SIZE=matrix%M*matrix%N
          matrix%NUMBER_NON_ZEROS=matrix%M*matrix%N
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=matrix%N
        CASE(matrix_diagonal_storage_type)
          IF(matrix%MAX_M==-1) matrix%MAX_M=matrix%M
          IF(matrix%MAX_N==-1) matrix%MAX_N=matrix%N
          matrix%SIZE=matrix%M
          matrix%NUMBER_NON_ZEROS=matrix%M
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=1
        CASE(matrix_column_major_storage_type)
          IF(matrix%MAX_M==-1) CALL flagerror("Maximum number of rows has not been set for this matrix.",err,error,*999)
          IF(matrix%MAX_N==-1) CALL flagerror("Maximum number of columns has not been set for this matrix.",err,error,*999)
          matrix%SIZE=matrix%MAX_M*matrix%N
          matrix%NUMBER_NON_ZEROS=matrix%M*matrix%N
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=matrix%N
        CASE(matrix_row_major_storage_type)
          IF(matrix%MAX_M==-1) CALL flagerror("Maximum number of rows has not been set for this matrix.",err,error,*999)
          IF(matrix%MAX_N==-1) CALL flagerror("Maximum number of columns has not been set for this matrix.",err,error,*999)
          matrix%SIZE=matrix%M*matrix%MAX_N
          matrix%NUMBER_NON_ZEROS=matrix%M*matrix%N
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=matrix%N
        CASE(matrix_compressed_row_storage_type)
          IF(matrix%NUMBER_NON_ZEROS==-1) CALL flagerror("Number of non-zeros has not been set for this matrix.",err,error,*999)
          IF(matrix%MAX_M==-1) matrix%MAX_M=matrix%M
          IF(matrix%MAX_N==-1) matrix%MAX_N=matrix%N
          matrix%SIZE=matrix%NUMBER_NON_ZEROS
          IF(.NOT.ALLOCATED(matrix%COLUMN_INDICES))  &
            & CALL flagerror("Matrix storage locations column indices have not been set.",err,error,*999)
          IF(.NOT.ALLOCATED(matrix%ROW_INDICES))  &
            & CALL flagerror("Matrix storage locations row indices have not been set.",err,error,*999)
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=0
          DO row_idx=1,matrix%M
            IF((matrix%ROW_INDICES(row_idx+1)-matrix%ROW_INDICES(row_idx))>matrix%MAXIMUM_COLUMN_INDICES_PER_ROW) &
              & matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=matrix%ROW_INDICES(row_idx+1)-matrix%ROW_INDICES(row_idx)
          ENDDO !row_idx
        CASE(matrix_compressed_column_storage_type)
          IF(matrix%NUMBER_NON_ZEROS==-1) CALL flagerror("Number of non-zeros has not been set for this matrix.",err,error,*999)
          IF(matrix%MAX_M==-1) matrix%MAX_M=matrix%M
          IF(matrix%MAX_N==-1) matrix%MAX_N=matrix%N
          matrix%SIZE=matrix%NUMBER_NON_ZEROS
          IF(.NOT.ALLOCATED(matrix%COLUMN_INDICES))  &
            & CALL flagerror("Matrix storage locations column indices have not been set.",err,error,*999)
          IF(.NOT.ALLOCATED(matrix%ROW_INDICES))  &
            & CALL flagerror("Matrix storage locations row indices have not been set.",err,error,*999)
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=0
          DO row_idx=1,matrix%M
            count=0
            DO column_idx=1,matrix%N
              DO row_idx2=matrix%COLUMN_INDICES(column_idx),matrix%COLUMN_INDICES(column_idx+1)-1
                IF(matrix%ROW_INDICES(row_idx2)==row_idx) count=count+1
              ENDDO !row_idx2
            ENDDO !column_idx
            IF(count>matrix%MAXIMUM_COLUMN_INDICES_PER_ROW) matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=count
          ENDDO !row_idx
        CASE(matrix_row_column_storage_type)
          IF(matrix%NUMBER_NON_ZEROS==-1) CALL flagerror("Number of non-zeros has not been set for this matrix.",err,error,*999)
          IF(matrix%MAX_M==-1) matrix%MAX_M=matrix%M
          IF(matrix%MAX_N==-1) matrix%MAX_N=matrix%N
          matrix%SIZE=matrix%NUMBER_NON_ZEROS  
          IF(.NOT.ALLOCATED(matrix%COLUMN_INDICES))  &
            & CALL flagerror("Matrix storage locations column indices have not been set.",err,error,*999)
          IF(.NOT.ALLOCATED(matrix%ROW_INDICES))  &
            & CALL flagerror("Matrix storage locations row indices have not been set.",err,error,*999)
          matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=0
          DO row_idx=1,matrix%M
            count=0
            DO row_idx2=1,matrix%NUMBER_NON_ZEROS
              IF(matrix%ROW_INDICES(row_idx2)==row_idx) count=count+1
            ENDDO !row_idx2            
            IF(count>matrix%MAXIMUM_COLUMN_INDICES_PER_ROW) matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=count
          ENDDO !row_idx          
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
        IF(matrix%SIZE>0) THEN
          SELECT CASE(matrix%DATA_TYPE)
          CASE(matrix_vector_intg_type)
            ALLOCATE(matrix%DATA_INTG(matrix%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate matrix integer data.",err,error,*999)
          CASE(matrix_vector_sp_type)
            ALLOCATE(matrix%DATA_SP(matrix%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate matrix single precision data.",err,error,*999)
          CASE(matrix_vector_dp_type)
            ALLOCATE(matrix%DATA_DP(matrix%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate matrix double precision data.",err,error,*999)
          CASE(matrix_vector_l_type)
            ALLOCATE(matrix%DATA_L(matrix%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate matrix logical data.",err,error,*999)
          CASE DEFAULT
            local_error="The matrix data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))//" is invalid."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        ENDIF
        matrix%ID=matrix_vector_id
        matrix_vector_id=matrix_vector_id+1
        matrix%MATRIX_FINISHED=.true.
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_CREATE_FINISH")
    RETURN
!!TODO: deallocate on error
999 errorsexits("MATRIX_CREATE_FINISH",err,error)
    RETURN 1
  END SUBROUTINE matrix_create_finish

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_create_start(MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("MATRIX_CREATE_START",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      CALL flagerror("Matrix is already associated.",err,error,*998)
    ELSE
      ALLOCATE(matrix,stat=err)
      IF(err/=0) CALL flagerror("Could not allocate the matrix.",err,error,*999)
      CALL matrix_initialise(matrix,err,error,*999)
      !Set the defaults
      matrix%DATA_TYPE=matrix_vector_dp_type
      matrix%STORAGE_TYPE=matrix_block_storage_type
    ENDIF
    
    exits("MATRIX_CREATE_START")
    RETURN
999 IF(ASSOCIATED(matrix)) CALL matrix_finalise(matrix,err,error,*998)
998 errorsexits("MATRIX_CREATE_START",err,error)
    RETURN 1
  END SUBROUTINE matrix_create_start

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_data_get_intg(MATRIX,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_DATA_GET_INTG",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)    
        IF(matrix%MATRIX_FINISHED) THEN
          IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
            data=>matrix%DATA_INTG
          ELSE
            local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the integer data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The matrix has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_DATA_GET_INTG")
    RETURN
999 errorsexits("MATRIX_DATA_GET_INTG",err,error)
    RETURN 1
  END SUBROUTINE matrix_data_get_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_data_get_sp(MATRIX,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    REAL(SP), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_DATA_GET_SP",err,error,*999)

     IF(ASSOCIATED(matrix)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(matrix%MATRIX_FINISHED) THEN
          IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
            data=>matrix%DATA_SP
          ELSE
            local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the single precision data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The matrix has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_DATA_GET_SP")
    RETURN
999 errorsexits("MATRIX_DATA_GET_SP",err,error)
    RETURN 1
  END SUBROUTINE matrix_data_get_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_data_get_dp(MATRIX,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    REAL(DP), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_DATA_GET_DP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(matrix%MATRIX_FINISHED) THEN
          IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
            data=>matrix%DATA_DP
          ELSE
            local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the double precision data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The matrix has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_DATA_GET_DP")
    RETURN
999 errorsexits("MATRIX_DATA_GET_DP",err,error)
    RETURN 1
  END SUBROUTINE matrix_data_get_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_data_get_l(MATRIX,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    LOGICAL, POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_DATA_GET_L",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(matrix%MATRIX_FINISHED) THEN
          IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
            data=>matrix%DATA_L
          ELSE
            local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the logical data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The matrix has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_DATA_GET_L")
    RETURN
999 errorsexits("MATRIX_DATA_GET_L",err,error)
    RETURN 1
  END SUBROUTINE matrix_data_get_l

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_datatypeget(matrix,dataType,err,error,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: matrix
    INTEGER(INTG), INTENT(OUT) :: dataType
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("Matrix_DataTypeGet",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(.NOT.matrix%matrix_finished) THEN
        CALL flag_error("The matrix has not been finished.",err,error,*999)
      ELSE
        datatype=matrix%data_type
      END IF
    ELSE
      CALL flag_error("Matrix is not associated.",err,error,*999)
    END IF

    exits("Matrix_DataTypeGet")
    RETURN
999 errorsexits("Matrix_DataTypeGet",err,error)
    RETURN 1
  END SUBROUTINE matrix_datatypeget

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_data_type_set(MATRIX,DATA_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: DATA_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_DATA_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN      
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has been finished.",err,error,*999)
      ELSE
        SELECT CASE(data_type)
        CASE(matrix_vector_intg_type)
          matrix%DATA_TYPE=matrix_vector_intg_type
        CASE(matrix_vector_sp_type)
          matrix%DATA_TYPE=matrix_vector_sp_type
        CASE(matrix_vector_dp_type)
          matrix%DATA_TYPE=matrix_vector_dp_type
        CASE(matrix_vector_l_type)
          matrix%DATA_TYPE=matrix_vector_l_type
        CASE DEFAULT
          local_error="The matrix vector data type of "//trim(number_to_vstring(data_type,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_DATA_TYPE_SET")
    RETURN
999 errorsexits("MATRIX_DATA_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_data_type_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_destroy(MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    enters("MATRIX_DESTROY",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      CALL matrix_finalise(matrix,err,error,*999)
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_DESTROY")
    RETURN
999 errorsexits("MATRIX_DESTROY",err,error)
    RETURN 1
  END SUBROUTINE matrix_destroy

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_duplicate(MATRIX,NEW_MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    TYPE(matrix_type), POINTER :: NEW_MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("MATRIX_DUPLICATE",err,error,*998)

    IF(ASSOCIATED(matrix)) THEN
      IF(ASSOCIATED(new_matrix)) THEN
        CALL flagerror("New matrix is already associated.",err,error,*998)
      ELSE
        CALL matrix_create_start(new_matrix,err,error,*999)
        CALL matrix_data_type_set(new_matrix,matrix%DATA_TYPE,err,error,*999)
        CALL matrix_size_set(new_matrix,matrix%M,matrix%N,err,error,*999)
        CALL matrix_storage_type_set(new_matrix,matrix%STORAGE_TYPE,err,error,*999)
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type,matrix_diagonal_storage_type)
          !Do nothing
        CASE(matrix_column_major_storage_type,matrix_row_major_storage_type)
          CALL matrix_max_size_set(new_matrix,matrix%MAX_M,matrix%MAX_N,err,error,*999)          
        CASE(matrix_compressed_row_storage_type,matrix_compressed_column_storage_type,matrix_row_column_storage_type)
          CALL matrix_number_non_zeros_set(new_matrix,matrix%NUMBER_NON_ZEROS,err,error,*999)
          CALL matrix_storage_locations_set(new_matrix,matrix%ROW_INDICES,matrix%COLUMN_INDICES,err,error,*999)
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
        CALL matrix_create_finish(new_matrix,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*998)
    ENDIF

    exits("MATRIX_DUPLICATE")
    RETURN
999 CALL matrix_finalise(new_matrix,err,error,*998)
998 errorsexits("MATRIX_DUPLICATE",err,error)
    RETURN 1
  END SUBROUTINE matrix_duplicate

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_finalise(MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("MATRIX_FINALISE",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(ALLOCATED(matrix%ROW_INDICES)) DEALLOCATE(matrix%ROW_INDICES)
      IF(ALLOCATED(matrix%COLUMN_INDICES)) DEALLOCATE(matrix%COLUMN_INDICES)
      IF(ALLOCATED(matrix%DATA_INTG)) DEALLOCATE(matrix%DATA_INTG)
      IF(ALLOCATED(matrix%DATA_SP)) DEALLOCATE(matrix%DATA_SP)
      IF(ALLOCATED(matrix%DATA_DP)) DEALLOCATE(matrix%DATA_DP)
      IF(ALLOCATED(matrix%DATA_L)) DEALLOCATE(matrix%DATA_L)
      DEALLOCATE(matrix)
    ENDIF

    exits("MATRIX_FINALISE")
    RETURN
999 errorsexits("MATRIX_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE matrix_finalise

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_initialise(MATRIX,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("MATRIX_INITIALISE",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      !!TODO: have a matrix user number etc.
      matrix%ID=0
      matrix%MATRIX_FINISHED=.false.
      matrix%M=0
      matrix%N=0
      matrix%MAX_M=-1
      matrix%MAX_N=-1
      matrix%DATA_TYPE=0
      matrix%STORAGE_TYPE=0
      matrix%NUMBER_NON_ZEROS=0
      matrix%SIZE=0      
      matrix%MAXIMUM_COLUMN_INDICES_PER_ROW=0      
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_INITIALISE")
    RETURN
999 errorsexits("MATRIX_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE matrix_initialise

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_max_columns_per_row_get(MATRIX,MAX_COLUMNS_PER_ROW,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: MAX_COLUMNS_PER_ROW
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("MATRIX_MAX_COLUMNS_PER_ROW_GET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        max_columns_per_row=matrix%MAXIMUM_COLUMN_INDICES_PER_ROW
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_MAX_COLUMNS_PER_ROW_GET")
    RETURN
999 errorsexits("MATRIX_MAX_COLUMNS_PER_ROW_GET",err,error)
    RETURN 1
  END SUBROUTINE matrix_max_columns_per_row_get

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_number_non_zeros_set(MATRIX,NUMBER_NON_ZEROS,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: NUMBER_NON_ZEROS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_NUMBER_NON_ZEROS_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has already been finished.",err,error,*999)
      ELSE
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          CALL flagerror("Can not set the number of non-zeros for a matrix with block storage.",err,error,*999)
        CASE(matrix_diagonal_storage_type)
          CALL flagerror("Can not set the number of non-zeros for a matrix with diagonal storage.",err,error,*999)
        CASE(matrix_column_major_storage_type)
          CALL flagerror("Can not set the number of non-zeros for a matrix with column major storage.",err,error,*999)          
        CASE(matrix_row_major_storage_type)
          CALL flagerror("Can not set the number of non-zeros for a matrix with row major storage.",err,error,*999)          
        CASE(matrix_compressed_row_storage_type,matrix_compressed_column_storage_type,matrix_row_column_storage_type)
          IF(number_non_zeros>=0) THEN
            matrix%NUMBER_NON_ZEROS=number_non_zeros
          ELSE
            local_error="The number of non-zeros ("//trim(number_to_vstring(number_non_zeros,"*",err,error))// &
              & ") is invalid. The number must be greater than or equal to zero."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_NUMBER_NON_ZEROS_SET")
    RETURN
999 errorsexits("MATRIX_NUMBER_NON_ZEROS_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_number_non_zeros_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_number_non_zeros_get(MATRIX,NUMBER_NON_ZEROS,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: NUMBER_NON_ZEROS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_NUMBER_NON_ZEROS_GET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type,matrix_diagonal_storage_type,matrix_column_major_storage_type, &
          & matrix_row_major_storage_type,matrix_compressed_row_storage_type,matrix_compressed_column_storage_type, &
          & matrix_row_column_storage_type)
          number_non_zeros=matrix%NUMBER_NON_ZEROS
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The matrix is not finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_NUMBER_NON_ZEROS_GET")
    RETURN
999 errorsexits("MATRIX_NUMBER_NON_ZEROS_GET",err,error)
    RETURN 1
  END SUBROUTINE matrix_number_non_zeros_get

  !
  !================================================================================================================================
  !
  SUBROUTINE matrix_linklist_set(MATRIX,LIST,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    TYPE(linkedlist),pointer :: LIST(:) 
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("MATRIX_LINKLIST_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has been finished",err,error,*999)
      ELSE
        matrix%LIST => list      
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_LINKLIST_SET")
    RETURN
999 errorsexits("MATRIX_LINKLIST_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_linklist_set

   !
  !================================================================================================================================
  !!>Gets the maximum number of columns in each row of a distributed matrix.
  SUBROUTINE matrix_linklist_get(MATRIX,LIST,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    type(linkedlist),pointer :: list(:) 
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    
    enters("MATRIX_LINKLIST_GET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        list=>matrix%LIST
      ELSE
        CALL flagerror("The matrix has not been finished",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_LINKLIST_GET")
    RETURN
999 errorsexits("MATRIX_LINKLIST_GET",err,error)
    RETURN 1
  END SUBROUTINE matrix_linklist_get

   !
  !================================================================================================================================
  !
  !
  SUBROUTINE matrix_max_size_set(MATRIX,MAX_M,MAX_N,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: MAX_M
    INTEGER(INTG), INTENT(IN) :: MAX_N
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_MAX_SIZE_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has been finished.",err,error,*999)
      ELSE
        IF(max_m>0) THEN
          IF(max_n>0) THEN
            IF(max_m>=matrix%M) THEN
              IF(max_n>=matrix%N) THEN
                matrix%MAX_M=max_m
                matrix%MAX_N=max_n
              ELSE
                local_error="The maximum number of matrix rows ("//trim(number_to_vstring(max_n,"*",err,error))// &
                  & ") must be >= the number of matrix rows ("//trim(number_to_vstring(matrix%N,"*",err,error))//")."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ELSE
              local_error="The maximum number of matrix columns ("//trim(number_to_vstring(max_m,"*",err,error))// &
                & ") must be >= the number of matrix columns ("//trim(number_to_vstring(matrix%M,"*",err,error))//")."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The maximum number of matrix columns of "//trim(number_to_vstring(max_n,"*",err,error))// &
              & " is invalid. The number must be > 0."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The maximum number of matrix rows of "//trim(number_to_vstring(max_m,"*",err,error))// &
            & " is invalid. The number must be > 0."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_MAX_SIZE_SET")
    RETURN
999 errorsexits("MATRIX_MAX_SIZE_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_max_size_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_output(ID,MATRIX,ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: ID
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: i,j
    CHARACTER(LEN=9) :: ROW_STRING,COL_STRING
    CHARACTER(LEN=39) :: INITIAL_STRING
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_OUTPUT",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          SELECT CASE(matrix%DATA_TYPE)
          CASE(matrix_vector_intg_type)
            CALL write_string_matrix(id,1,1,matrix%M,1,1,matrix%N,8,8,reshape(matrix%DATA_INTG,(/matrix%MAX_M,matrix%MAX_N/)), &
              & write_string_matrix_name_and_indices,'("Matrix','(",I9,",:)',':",8(X,I13))','(20X,8(X,I13))', &
              & err,error,*999)
          CASE(matrix_vector_sp_type)
            CALL write_string_matrix(id,1,1,matrix%M,1,1,matrix%N,8,8,reshape(matrix%DATA_SP,(/matrix%MAX_M,matrix%MAX_N/)), &
              & write_string_matrix_name_and_indices,'("Matrix','(",I9,",:)',':",8(X,E13.6))','(20X,8(X,E13.6))', &
              & err,error,*999)
          CASE(matrix_vector_dp_type)
            CALL write_string_matrix(id,1,1,matrix%M,1,1,matrix%N,8,8,reshape(matrix%DATA_DP,(/matrix%MAX_M,matrix%MAX_N/)), &
              & write_string_matrix_name_and_indices,'("Matrix','(",I9,",:)',':",8(X,E13.6))','(20X,8(X,E13.6))', &
              & err,error,*999)
          CASE(matrix_vector_l_type)            
            CALL write_string_matrix(id,1,1,matrix%M,1,1,matrix%N,8,8,reshape(matrix%DATA_L,(/matrix%MAX_M,matrix%MAX_N/)), &
              & write_string_matrix_name_and_indices,'("Matrix','(",I9,",:)',':",8(X,L13))','(20X,8(X,L13))', &
              & err,error,*999)
          CASE DEFAULT
            local_error="The matrix data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))//" is invalid."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(matrix_diagonal_storage_type)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(matrix_column_major_storage_type)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(matrix_row_major_storage_type)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(matrix_compressed_row_storage_type)
          DO i=1,matrix%M
            row_string=number_to_character(i,"I9",err,error)
            SELECT CASE(matrix%DATA_TYPE)              
            CASE(matrix_vector_intg_type)
              initial_string='("Matrix('//row_string//',:):",8(X,I13))'
              CALL write_string_vector(id,matrix%ROW_INDICES(i),1,matrix%ROW_INDICES(i+1)-1,8,8,matrix%DATA_INTG,initial_string, &
                & '(20X,8(X,I13))',err,error,*999)
            CASE(matrix_vector_sp_type)
              initial_string='("Matrix('//row_string//',:):",8(X,E13.6))'
              CALL write_string_vector(id,matrix%ROW_INDICES(i),1,matrix%ROW_INDICES(i+1)-1,8,8,matrix%DATA_SP,initial_string, &
                & '(20X,8(X,E13.6))',err,error,*999)
            CASE(matrix_vector_dp_type)
              initial_string='("Matrix('//row_string//',:):",8(X,E13.6))'
              CALL write_string_vector(id,matrix%ROW_INDICES(i),1,matrix%ROW_INDICES(i+1)-1,8,8,matrix%DATA_DP,initial_string, &
                & '(20X,8(X,E13.6))',err,error,*999)
            CASE(matrix_vector_l_type)            
              initial_string='("Matrix('//row_string//',:):",8(X,L13))'
              CALL write_string_vector(id,matrix%ROW_INDICES(i),1,matrix%ROW_INDICES(i+1)-1,8,8,matrix%DATA_L,initial_string, &
                & '(20X,8(X,L13))',err,error,*999)
            CASE DEFAULT
              local_error="The matrix data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(local_error,err,error,*999)
            END SELECT
          ENDDO !i
        CASE(matrix_compressed_column_storage_type)
          DO j=1,matrix%N
            col_string=number_to_character(j,"I9",err,error)
            SELECT CASE(matrix%DATA_TYPE)              
            CASE(matrix_vector_intg_type)
              initial_string='("Matrix(:,'//col_string//'):",8(X,I13))'
              CALL write_string_vector(id,matrix%COLUMN_INDICES(j),1,matrix%COLUMN_INDICES(j+1)-1,8,8,matrix%DATA_INTG, &
                & initial_string,'(20X,8(X,I13))',err,error,*999)
            CASE(matrix_vector_sp_type)
              initial_string='("Matrix(:,'//col_string//'):",8(X,E13.6))'
              CALL write_string_vector(id,matrix%COLUMN_INDICES(j),1,matrix%COLUMN_INDICES(j+1)-1,8,8,matrix%DATA_SP, &
                & initial_string,'(20X,8(X,E13.6))',err,error,*999)
            CASE(matrix_vector_dp_type)
              initial_string='("Matrix(:,'//col_string//'):",8(X,E13.6))'
              CALL write_string_vector(id,matrix%COLUMN_INDICES(j),1,matrix%COLUMN_INDICES(j+1)-1,8,8,matrix%DATA_DP, &
                & initial_string,'(20X,8(X,E13.6))',err,error,*999)
            CASE(matrix_vector_l_type)            
              initial_string='("Matrix(:,'//col_string//'):",8(X,L13))'
              CALL write_string_vector(id,matrix%COLUMN_INDICES(j),1,matrix%COLUMN_INDICES(j+1)-1,8,8,matrix%DATA_L, &
                & initial_string,'(20X,8(X,L13))',err,error,*999)
            CASE DEFAULT
              local_error="The matrix data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))//" is invalid."
              CALL flagerror(local_error,err,error,*999)
            END SELECT
          ENDDO !j
        CASE(matrix_row_column_storage_type)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF      
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_OUTPUT")
    RETURN
999 errorsexits("MATRIX_OUTPUT",err,error)
    RETURN 1
  END SUBROUTINE matrix_output

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_size_set(MATRIX,M,N,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: M
    INTEGER(INTG), INTENT(IN) :: N
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_SIZE_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has been finished.",err,error,*999)
      ELSE
        IF(m>0) THEN
          IF(n>0) THEN
            matrix%M=m
            matrix%N=n
          ELSE
            local_error="The number of matrix columns of "//trim(number_to_vstring(n,"*",err,error))// &
              & " is invalid. The number must be >0."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The number of matrix rows of "//trim(number_to_vstring(m,"*",err,error))// &
            & " is invalid. The number must be >0."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_SIZE_SET")
    RETURN
999 errorsexits("MATRIX_SIZE_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_size_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_storage_location_find(MATRIX,I,J,LOCATION,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: I
    INTEGER(INTG), INTENT(IN) :: J
    INTEGER(INTG), INTENT(OUT) :: LOCATION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOWLIMIT,MIDPOINT,UPLIMIT
    LOGICAL :: FOUNDCOLUMN, FOUNDROW
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_STORAGE_LOCATION_FIND",err,error,*999)

    location=0
    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(i<1.OR.i>matrix%M) THEN
          local_error="Row number "//trim(number_to_vstring(i,"*",err,error))//" is outside the matrix range of 1 to "// &
            & trim(number_to_vstring(matrix%M,"*",err,error))//"."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
        IF(j<1.OR.j>matrix%N) THEN
          local_error="Column number "//trim(number_to_vstring(j,"*",err,error))//" is outside the matrix range of 1 to "// &
            & trim(number_to_vstring(matrix%M,"*",err,error))//"."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          location=i+(j-1)*matrix%M
        CASE(matrix_diagonal_storage_type)
          IF(i==j) location=i
        CASE(matrix_column_major_storage_type)
          location=i+(j-1)*matrix%MAX_M
        CASE(matrix_row_major_storage_type)
          location=(i-1)*matrix%MAX_N+j
        CASE(matrix_compressed_row_storage_type)
          !Search for the column number in the sparsity list using the bisection (binary search) algorithm
          lowlimit=matrix%ROW_INDICES(i)
          IF(j>=matrix%COLUMN_INDICES(lowlimit)) THEN
            uplimit=matrix%ROW_INDICES(i+1)
            IF(uplimit>lowlimit) THEN
              IF(j<=matrix%COLUMN_INDICES(uplimit-1)) THEN
                DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                  midpoint=(uplimit+lowlimit)/2
                  IF(matrix%COLUMN_INDICES(midpoint)>j) THEN
                    uplimit=midpoint
                  ELSE
                    lowlimit=midpoint
                  ENDIF
                ENDDO
                DO k=lowlimit,uplimit
                  IF(matrix%COLUMN_INDICES(k)==j) THEN
                    location=k
                    EXIT
                  ENDIF
                ENDDO !k
              ENDIF
            ENDIF
          ENDIF
        CASE(matrix_compressed_column_storage_type)
          !Search for the row number in the sparsity list using the bisection (binary search) algorithm
          lowlimit=matrix%COLUMN_INDICES(j)
          IF(i>=matrix%ROW_INDICES(lowlimit)) THEN
            uplimit=matrix%COLUMN_INDICES(j+1)
            IF(uplimit>lowlimit) THEN
              IF(i<=matrix%ROW_INDICES(uplimit-1)) THEN
                DO WHILE((uplimit-lowlimit)>1)
                  midpoint=(uplimit+lowlimit)/2
                  IF(matrix%ROW_INDICES(midpoint)>i) THEN
                    uplimit=midpoint
                  ELSE
                    lowlimit=midpoint
                  ENDIF
                ENDDO
                IF(matrix%ROW_INDICES(lowlimit)==i) location=lowlimit
                DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                  midpoint=(uplimit+lowlimit)/2
                  IF(matrix%ROW_INDICES(midpoint)>i) THEN
                    uplimit=midpoint
                  ELSE
                    lowlimit=midpoint
                  ENDIF
                ENDDO
                DO k=lowlimit,uplimit
                  IF(matrix%ROW_INDICES(k)==i) THEN
                    location=k
                    EXIT
                  ENDIF
                ENDDO !k
              ENDIF
            ENDIF
          ENDIF
        CASE(matrix_row_column_storage_type)
          foundrow=.false.
          location=1
          DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
            IF(matrix%ROW_INDICES(location)==i) THEN
              DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                IF(matrix%COLUMN_INDICES(location)==j.AND.matrix%ROW_INDICES(location)==i) THEN
                  foundcolumn=.true.
                ELSE IF(matrix%ROW_INDICES(location)/=i) THEN
                  location=matrix%SIZE+1
                ELSE
                  location=location+1
                ENDIF
              ENDDO
            ELSE
              location=location+1
            ENDIF
          ENDDO
          IF(.NOT.(foundrow.AND.foundcolumn)) location=0
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_STORAGE_LOCATION_FIND")
    RETURN
999 errorsexits("MATRIX_STORAGE_LOCATION_FIND",err,error)
    RETURN 1
  END SUBROUTINE matrix_storage_location_find

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_storage_locations_get(MATRIX,ROW_INDICES,COLUMN_INDICES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), POINTER :: ROW_INDICES(:)
    INTEGER(INTG), POINTER :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_STORAGE_LOCATIONS_GET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          CALL flagerror("Can not get matrix locations for a block storage matrix.",err,error,*999)
        CASE(matrix_diagonal_storage_type)
          CALL flagerror("Can not get matrix locations for a diagonal storage matrix.",err,error,*999)
        CASE(matrix_column_major_storage_type)
          CALL flagerror("Can not get matrix locations for a column major storage matrix.",err,error,*999)
        CASE(matrix_row_major_storage_type)
          CALL flagerror("Can not get matrix locations for a row major storage matrix.",err,error,*999)
        CASE(matrix_compressed_row_storage_type)          
          row_indices=>matrix%ROW_INDICES
          column_indices=>matrix%COLUMN_INDICES
        CASE(matrix_compressed_column_storage_type)
          row_indices=>matrix%ROW_INDICES
          column_indices=>matrix%COLUMN_INDICES          
        CASE(matrix_row_column_storage_type)
          row_indices=>matrix%ROW_INDICES
          column_indices=>matrix%COLUMN_INDICES
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Matrix has not been finished.",err,error,*999)
      ENDIF
   ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_STORAGE_LOCATIONS_GET")
    RETURN
999 errorsexits("MATRIX_STORAGE_LOCATIONS_GET",err,error)
    RETURN 1
  END SUBROUTINE matrix_storage_locations_get

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_storage_locations_set(MATRIX,ROW_INDICES,COLUMN_INDICES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_STORAGE_LOCATIONS_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("Matrix has been finished.",err,error,*999)
      ELSE
        SELECT CASE(matrix%STORAGE_TYPE)
        CASE(matrix_block_storage_type)
          CALL flagerror("Can not set matrix locations for a block storage matrix.",err,error,*999)
        CASE(matrix_diagonal_storage_type)
          CALL flagerror("Can not set matrix locations for a diagonal storage matrix.",err,error,*999)
        CASE(matrix_column_major_storage_type)
          CALL flagerror("Can not set matrix locations for a column major storage matrix.",err,error,*999)
        CASE(matrix_row_major_storage_type)
          CALL flagerror("Can not set matrix locations for a row major storage matrix.",err,error,*999)
        CASE(matrix_compressed_row_storage_type)
          IF(SIZE(row_indices,1)==matrix%M+1) THEN
            IF(SIZE(column_indices,1)==matrix%NUMBER_NON_ZEROS) THEN
              IF(row_indices(1)==1) THEN
                IF(row_indices(matrix%M+1)==matrix%NUMBER_NON_ZEROS+1) THEN
                  DO i=2,matrix%M+1
                    IF(row_indices(i)<row_indices(i-1)) THEN
                      local_error="Invalid row indices. Row "//trim(number_to_vstring(i,"*",err,error))//" index number ("// &
                        & trim(number_to_vstring(row_indices(i),"*",err,error))//") is less than row "// &
                        & trim(number_to_vstring(i-1,"*",err,error))//" index number ("// &
                        & trim(number_to_vstring(row_indices(i-1),"*",err,error))//")."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF                    
                  ENDDO !i
                  DO i=1,matrix%M
                    DO j=row_indices(i),row_indices(i+1)-1
                      k=column_indices(j)
                      IF(k>0) THEN
                        IF(k>matrix%N) THEN
                          local_error="Invalid column indices. Column index "//trim(number_to_vstring(j,"*",err,error))//" ("// &
                            & trim(number_to_vstring(k,"*",err,error))//") is greater than the number of columns ("// &
                            & trim(number_to_vstring(matrix%N,"*",err,error))//")."
                          CALL flagerror(local_error,err,error,*999)
                        ENDIF
                      ELSE
                        local_error="Invalid column indices. Column index "//trim(number_to_vstring(j,"*",err,error))//" ("// &
                          & trim(number_to_vstring(k,"*",err,error))//") is less than zero."
                        CALL flagerror(local_error,err,error,*999)
                      ENDIF
                    ENDDO !j
                  ENDDO !i
                  IF(ALLOCATED(matrix%ROW_INDICES)) DEALLOCATE(matrix%ROW_INDICES)
                  IF(ALLOCATED(matrix%COLUMN_INDICES)) DEALLOCATE(matrix%COLUMN_INDICES)
                  ALLOCATE(matrix%ROW_INDICES(matrix%M+1),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate matrix row indices.",err,error,*999)
                  ALLOCATE(matrix%COLUMN_INDICES(matrix%NUMBER_NON_ZEROS),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate matrix column indices.",err,error,*999)                  
                  matrix%ROW_INDICES(1:matrix%M+1)=row_indices(1:matrix%M+1)
                  matrix%COLUMN_INDICES(1:matrix%NUMBER_NON_ZEROS)=column_indices(1:matrix%NUMBER_NON_ZEROS)
                  !Don't really need this???
                  !DO i=1,MATRIX%M
                  !  CALL LIST_SORT(MATRIX%COLUMN_INDICES(MATRIX%ROW_INDICES(i):MATRIX%ROW_INDICES(i+1)-1),ERR,ERROR,*999)
                  !ENDDO !i
                ELSE
                  local_error="Invalid row indices. The last row index ("// &
                    & trim(number_to_vstring(row_indices(matrix%M+1),"*",err,error))// &
                    & ") does not equal the number of non-zeros + 1 ("// &
                    & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS+1,"*",err,error))//")."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ELSE
                local_error="Invalid row indices. The first row index ("// &
                  & trim(number_to_vstring(row_indices(1),"*",err,error))//") does not equal 1."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ELSE
              local_error="The supplied number of column indices ("// &
                & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
                & ") does not match the number of non-zeros in the matrix ("// &
                & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS,"*",err,error))//")."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The supplied number of row indices ("//trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
              & ") does not match the number of rows in the matrix + 1 ("// &
              & trim(number_to_vstring(matrix%M+1,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        CASE(matrix_compressed_column_storage_type)
          IF(SIZE(column_indices,1)==matrix%N+1) THEN
            IF(SIZE(row_indices,1)==matrix%NUMBER_NON_ZEROS) THEN
             IF(column_indices(1)==1) THEN
                IF(column_indices(matrix%N+1)==matrix%NUMBER_NON_ZEROS+1) THEN
                  IF(column_indices(1)/=1) THEN
                    local_error="Invalid column indices. Column index 1 ("// &
                      & trim(number_to_vstring(column_indices(1),"*",err,error))//") "// &
                      & " should be equal to one."
                    CALL flagerror(local_error,err,error,*999)
                  END IF
                  DO j=2,matrix%N+1
                    IF(column_indices(j)<column_indices(j-1)) THEN
                      local_error="Invalid column indices. Column "//trim(number_to_vstring(j,"*",err,error))// &
                        & " index number ("//trim(number_to_vstring(column_indices(j),"*",err,error))//") is less than column "// &
                        & trim(number_to_vstring(j-1,"*",err,error))//" index number ("// &
                        & trim(number_to_vstring(column_indices(j-1),"*",err,error))//")."
                      CALL flagerror(local_error,err,error,*999)
                    END IF
                    IF(column_indices(j)<0.OR.column_indices(j)>matrix%NUMBER_NON_ZEROS+1) THEN
                      local_error="Invalid column indices. Column index "//trim(number_to_vstring(j,"*",err,error))//" ("// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//") "// &
                        & " should be in the range of one to the number of non-zeros + 1 ("// &
                        & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS+1,"*",err,error))//")."
                      CALL flagerror(local_error,err,error,*999)
                    END IF
                  ENDDO !i
                  DO j=1,matrix%N
                    DO i=column_indices(j),column_indices(j+1)-1
                      k=row_indices(i)
                      IF(k>0) THEN
                        IF(k>matrix%M) THEN
                          local_error="Invalid row indices. Row index "//trim(number_to_vstring(i,"*",err,error))//" ("// &
                            & trim(number_to_vstring(k,"*",err,error))//") is greater than the number of rows ("// &
                            & trim(number_to_vstring(matrix%M,"*",err,error))//")."
                          CALL flagerror(local_error,err,error,*999)
                        ENDIF
                      ELSE
                        local_error="Invalid row indices. Row index "//trim(number_to_vstring(i,"*",err,error))//" ("// &
                          & trim(number_to_vstring(k,"*",err,error))//") is less than zero."
                        CALL flagerror(local_error,err,error,*999)
                      ENDIF
                    ENDDO !i
                  ENDDO !j
                  IF(ALLOCATED(matrix%ROW_INDICES)) DEALLOCATE(matrix%ROW_INDICES)
                  IF(ALLOCATED(matrix%COLUMN_INDICES)) DEALLOCATE(matrix%COLUMN_INDICES)
                  ALLOCATE(matrix%ROW_INDICES(matrix%NUMBER_NON_ZEROS),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate matrix row indices.",err,error,*999)
                  ALLOCATE(matrix%COLUMN_INDICES(matrix%N+1),stat=err)
                  IF(err/=0) CALL flagerror("Could not allocate matrix column indices.",err,error,*999)                  
                  matrix%ROW_INDICES(1:matrix%NUMBER_NON_ZEROS)=row_indices(1:matrix%NUMBER_NON_ZEROS)
                  matrix%COLUMN_INDICES(1:matrix%N+1)=column_indices(1:matrix%N+1)
                  !Don't really need this???
                  !DO j=1,MATRIX%N                    
                  !  CALL LIST_SORT(MATRIX%ROW_INDICES(MATRIX%COLUMN_INDICES(j):MATRIX%COLUMN_INDICES(j+1)-1),ERR,ERROR,*999)
                  !ENDDO !j
                ELSE
                  local_error="Invalid column indices. The last column index ("// &
                    & trim(number_to_vstring(column_indices(matrix%N+1),"*",err,error))// &
                    & ") does not equal the number of non-zeros + 1 ("// &
                    & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS+1,"*",err,error))//")."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ELSE
                local_error="Invalid column indices. The first column index ("// &
                  & trim(number_to_vstring(column_indices(1),"*",err,error))//") does not equal 1."
                CALL flagerror(local_error,err,error,*999)
              ENDIF
            ELSE
              local_error="The supplied number of row indices ("// &
                & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
                & ") does not match the number of non-zeros in the matrix ("// &
                & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS,"*",err,error))//")."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The supplied number of column indices ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not match the number of columns in the matrix + 1 ("// &
              & trim(number_to_vstring(matrix%N+1,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        CASE(matrix_row_column_storage_type)
          IF(SIZE(row_indices,1)==matrix%NUMBER_NON_ZEROS) THEN
            IF(SIZE(column_indices,1)==matrix%NUMBER_NON_ZEROS) THEN
              DO k=1,matrix%NUMBER_NON_ZEROS
                IF(row_indices(k)<1.OR.row_indices(k)>matrix%M) THEN
                  local_error="Invalid row indices. Row index number "//trim(number_to_vstring(k,"*",err,error))//" ("// &
                    & trim(number_to_vstring(row_indices(k),"*",err,error))// &
                    & ") is out of range. The row index must be between 1 and "// &
                    & trim(number_to_vstring(matrix%M,"*",err,error))//"."
                  CALL flagerror(local_error,err,error,*999)
                ELSE IF(column_indices(k)<1.OR.column_indices(k)>matrix%N) THEN
                  local_error="Invalid column indices. Column index number "//trim(number_to_vstring(k,"*",err,error))//" ("// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))// &
                    & ") is out of range. The column index must be between 1 and "// &
                    & trim(number_to_vstring(matrix%N,"*",err,error))//"."
                  CALL flagerror(local_error,err,error,*999)
                ENDIF
              ENDDO !k
              matrix%ROW_INDICES(1:matrix%NUMBER_NON_ZEROS)=row_indices(1:matrix%NUMBER_NON_ZEROS)
              matrix%COLUMN_INDICES(1:matrix%NUMBER_NON_ZEROS)=column_indices(1:matrix%NUMBER_NON_ZEROS)
              !!TODO: sort the row and colum indices!!!!!
            ELSE
              local_error="The supplied number of column indices ("// &
                & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
                & ") does not match the number of non-zeros in the matrix ("// &
                & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS,"*",err,error))//")."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The supplied number of row indices ("// &
              & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
              & ") does not match the number of non-zeros in the matrix ("// &
              & trim(number_to_vstring(matrix%NUMBER_NON_ZEROS,"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_STORAGE_LOCATIONS_SET")
    RETURN
999 errorsexits("MATRIX_STORAGE_LOCATIONS_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_storage_locations_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_storage_type_get(MATRIX,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(OUT) :: STORAGE_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("MATRIX_STORAGE_TYPE_GET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN        
        storage_type=matrix%STORAGE_TYPE
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_STORAGE_TYPE_GET")
    RETURN
999 errorsexits("MATRIX_STORAGE_TYPE_GET",err,error)
    RETURN 1
  END SUBROUTINE matrix_storage_type_get

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_storage_type_set(MATRIX,STORAGE_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: STORAGE_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_STORAGE_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        CALL flagerror("The matrix has been finished.",err,error,*999)
      ELSE
        SELECT CASE(storage_type)
        CASE(matrix_block_storage_type)
          matrix%STORAGE_TYPE=matrix_block_storage_type
        CASE(matrix_diagonal_storage_type)
          matrix%STORAGE_TYPE=matrix_diagonal_storage_type
        CASE(matrix_column_major_storage_type)
          matrix%STORAGE_TYPE=matrix_column_major_storage_type
        CASE(matrix_row_major_storage_type)
          matrix%STORAGE_TYPE=matrix_row_major_storage_type
        CASE(matrix_compressed_row_storage_type)
          matrix%STORAGE_TYPE=matrix_compressed_row_storage_type
        CASE(matrix_compressed_column_storage_type)
          matrix%STORAGE_TYPE=matrix_compressed_column_storage_type
        CASE(matrix_row_column_storage_type)
          matrix%STORAGE_TYPE=matrix_row_column_storage_type
        CASE DEFAULT
          local_error="The matrix storage type of "//trim(number_to_vstring(storage_type,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_STORAGE_TYPE_SET")
    RETURN
999 errorsexits("MATRIX_STORAGE_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE matrix_storage_type_set

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_intg(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_INTG",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_INTG")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_INTG",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_intg1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    INTEGER(INTG), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_INTG1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_INTG1")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_INTG1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_intg1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_intg2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,k,ROW_INDEX,PREVIOUS_ROW_INDEX,COLUMN_INDEX,PREVIOUS_COLUMN_INDEX,LOCATION,LOWLIMIT,MIDPOINT,UPLIMIT
    LOGICAL :: FOUNDCOLUMN, FOUNDROW
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_INTG2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              SELECT CASE(matrix%STORAGE_TYPE)
              CASE(matrix_block_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_diagonal_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(row_index==column_index) location=row_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_column_major_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%MAX_M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_major_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    location=(row_index-1)*matrix%MAX_N+column_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_row_storage_type)
                !Search for the column number in the sparsity list using the bisection (binary search) algorithm
                previous_column_index=-1
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  lowlimit=matrix%ROW_INDICES(row_index)
                  uplimit=matrix%ROW_INDICES(row_index+1)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(column_index<=previous_column_index) THEN
                      lowlimit=matrix%ROW_INDICES(row_index)
                    ELSE
                      uplimit=matrix%ROW_INDICES(row_index+1)
                    ENDIF
                    previous_column_index=column_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%COLUMN_INDICES(midpoint)>column_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%COLUMN_INDICES(k)==column_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_column_storage_type)
                !Search for the row number in the sparsity list using the bisection (binary search) algorithm
                previous_row_index=-1
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  lowlimit=matrix%COLUMN_INDICES(column_index)
                  uplimit=matrix%COLUMN_INDICES(column_index+1)
                  DO i=1,SIZE(row_indices,1)
                    location=0
                    row_index=row_indices(i)
                    IF(row_index<=previous_row_index) THEN
                      lowlimit=matrix%COLUMN_INDICES(column_index)
                    ELSE
                      uplimit=matrix%COLUMN_INDICES(column_index+1)
                    ENDIF
                    previous_row_index=row_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%ROW_INDICES(midpoint)>row_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%ROW_INDICES(k)==row_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_column_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    foundrow=.false.
                    location=1
                    DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                      IF(matrix%ROW_INDICES(location)==row_index) THEN
                        DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                          IF(matrix%COLUMN_INDICES(location)==column_index.AND.matrix%ROW_INDICES(location)==row_index) THEN
                            foundcolumn=.true.
                          ELSE IF(matrix%ROW_INDICES(location)/=row_index) THEN
                            location=matrix%SIZE+1
                          ELSE
                            location=location+1
                          ENDIF
                        ENDDO
                      ELSE
                        location=location+1
                      ENDIF
                    ENDDO
                    IF(.NOT.(foundrow.AND.foundcolumn)) location=0
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_INTG(location)=matrix%DATA_INTG(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE DEFAULT
                local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_INTG2")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_INTG2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_intg2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_sp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_SP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_SP")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_SP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_sp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(SP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_SP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_SP(location)=matrix%DATA_SP(location)+VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_SP1")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_SP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_sp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_sp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,k,ROW_INDEX,PREVIOUS_ROW_INDEX,COLUMN_INDEX,PREVIOUS_COLUMN_INDEX,LOCATION,LOWLIMIT,MIDPOINT,UPLIMIT
    LOGICAL :: FOUNDCOLUMN, FOUNDROW
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_SP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              SELECT CASE(matrix%STORAGE_TYPE)
              CASE(matrix_block_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_diagonal_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(row_index==column_index) location=row_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_column_major_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%MAX_M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_major_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    location=(row_index-1)*matrix%MAX_N+column_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_row_storage_type)
                !Search for the column number in the sparsity list using the bisection (binary search) algorithm
                previous_column_index=-1
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  lowlimit=matrix%ROW_INDICES(row_index)
                  uplimit=matrix%ROW_INDICES(row_index+1)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(column_index<=previous_column_index) THEN
                      lowlimit=matrix%ROW_INDICES(row_index)
                    ELSE
                      uplimit=matrix%ROW_INDICES(row_index+1)
                    ENDIF
                    previous_column_index=column_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%COLUMN_INDICES(midpoint)>column_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%COLUMN_INDICES(k)==column_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_column_storage_type)
                !Search for the row number in the sparsity list using the bisection (binary search) algorithm
                previous_row_index=-1
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  lowlimit=matrix%COLUMN_INDICES(column_index)
                  uplimit=matrix%COLUMN_INDICES(column_index+1)
                  DO i=1,SIZE(row_indices,1)
                    location=0
                    row_index=row_indices(i)
                    IF(row_index<=previous_row_index) THEN
                      lowlimit=matrix%COLUMN_INDICES(column_index)
                    ELSE
                      uplimit=matrix%COLUMN_INDICES(column_index+1)
                    ENDIF
                    previous_row_index=row_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%ROW_INDICES(midpoint)>row_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%ROW_INDICES(k)==row_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_column_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    foundrow=.false.
                    location=1
                    DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                      IF(matrix%ROW_INDICES(location)==row_index) THEN
                        DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                          IF(matrix%COLUMN_INDICES(location)==column_index.AND.matrix%ROW_INDICES(location)==row_index) THEN
                            foundcolumn=.true.
                          ELSE IF(matrix%ROW_INDICES(location)/=row_index) THEN
                            location=matrix%SIZE+1
                          ELSE
                            location=location+1
                          ENDIF
                        ENDDO
                      ELSE
                        location=location+1
                      ENDIF
                    ENDDO
                    IF(.NOT.(foundrow.AND.foundcolumn)) location=0
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_SP(location)=matrix%DATA_SP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE DEFAULT
                local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_SP2")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_SP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_sp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_dp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_DP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_DP")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_DP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_dp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_DP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_DP(location)=matrix%DATA_DP(location)+VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_DP1")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_DP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_dp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_dp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,k,ROW_INDEX,PREVIOUS_ROW_INDEX,COLUMN_INDEX,PREVIOUS_COLUMN_INDEX,LOCATION,LOWLIMIT,MIDPOINT,UPLIMIT
    LOGICAL :: FOUNDCOLUMN, FOUNDROW
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_DP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              SELECT CASE(matrix%STORAGE_TYPE)
              CASE(matrix_block_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_diagonal_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(row_index==column_index) location=row_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_column_major_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%MAX_M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_major_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    location=(row_index-1)*matrix%MAX_N+column_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_row_storage_type)
                !Search for the column number in the sparsity list using the bisection (binary search) algorithm
                previous_column_index=-1
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  lowlimit=matrix%ROW_INDICES(row_index)
                  uplimit=matrix%ROW_INDICES(row_index+1)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(column_index<=previous_column_index) THEN
                      lowlimit=matrix%ROW_INDICES(row_index)
                    ELSE
                      uplimit=matrix%ROW_INDICES(row_index+1)
                    ENDIF
                    previous_column_index=column_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%COLUMN_INDICES(midpoint)>column_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%COLUMN_INDICES(k)==column_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_column_storage_type)
                !Search for the row number in the sparsity list using the bisection (binary search) algorithm
                previous_row_index=-1
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  lowlimit=matrix%COLUMN_INDICES(column_index)
                  uplimit=matrix%COLUMN_INDICES(column_index+1)
                  DO i=1,SIZE(row_indices,1)
                    location=0
                    row_index=row_indices(i)
                    IF(row_index<=previous_row_index) THEN
                      lowlimit=matrix%COLUMN_INDICES(column_index)
                    ELSE
                      uplimit=matrix%COLUMN_INDICES(column_index+1)
                    ENDIF
                    previous_row_index=row_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%ROW_INDICES(midpoint)>row_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%ROW_INDICES(k)==row_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_column_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    foundrow=.false.
                    location=1
                    DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                      IF(matrix%ROW_INDICES(location)==row_index) THEN
                        DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                          IF(matrix%COLUMN_INDICES(location)==column_index.AND.matrix%ROW_INDICES(location)==row_index) THEN
                            foundcolumn=.true.
                          ELSE IF(matrix%ROW_INDICES(location)/=row_index) THEN
                            location=matrix%SIZE+1
                          ELSE
                            location=location+1
                          ENDIF
                        ENDDO
                      ELSE
                        location=location+1
                      ENDIF
                    ENDDO
                    IF(.NOT.(foundrow.AND.foundcolumn)) location=0
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_DP(location)=matrix%DATA_DP(location)+values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE DEFAULT
                local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_DP2")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_DP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_dp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_l(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_L",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_L")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_L",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_l

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_l1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    LOGICAL, INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_L1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_L(location)=matrix%DATA_L(location).OR.VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_L1")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_L1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_l1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_add_l2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,k,ROW_INDEX,PREVIOUS_ROW_INDEX,COLUMN_INDEX,PREVIOUS_COLUMN_INDEX,LOCATION,LOWLIMIT,MIDPOINT,UPLIMIT
    LOGICAL :: FOUNDCOLUMN, FOUNDROW
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_ADD_L2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              SELECT CASE(matrix%STORAGE_TYPE)
              CASE(matrix_block_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_diagonal_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(row_index==column_index) location=row_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_column_major_storage_type)
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  DO i=1,SIZE(row_indices,1)
                    row_index=row_indices(i)
                    location=row_index+(column_index-1)*matrix%MAX_M
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_major_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    location=(row_index-1)*matrix%MAX_N+column_index
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_row_storage_type)
                !Search for the column number in the sparsity list using the bisection (binary search) algorithm
                previous_column_index=-1
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  lowlimit=matrix%ROW_INDICES(row_index)
                  uplimit=matrix%ROW_INDICES(row_index+1)
                  DO j=1,SIZE(column_indices,1)
                    location=0
                    column_index=column_indices(j)
                    IF(column_index<=previous_column_index) THEN
                      lowlimit=matrix%ROW_INDICES(row_index)
                    ELSE
                      uplimit=matrix%ROW_INDICES(row_index+1)
                    ENDIF
                    previous_column_index=column_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%COLUMN_INDICES(midpoint)>column_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%COLUMN_INDICES(k)==column_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE(matrix_compressed_column_storage_type)
                !Search for the row number in the sparsity list using the bisection (binary search) algorithm
                previous_row_index=-1
                DO j=1,SIZE(column_indices,1)
                  column_index=column_indices(j)
                  lowlimit=matrix%COLUMN_INDICES(column_index)
                  uplimit=matrix%COLUMN_INDICES(column_index+1)
                  DO i=1,SIZE(row_indices,1)
                    location=0
                    row_index=row_indices(i)
                    IF(row_index<=previous_row_index) THEN
                      lowlimit=matrix%COLUMN_INDICES(column_index)
                    ELSE
                      uplimit=matrix%COLUMN_INDICES(column_index+1)
                    ENDIF
                    previous_row_index=row_index
                    DO WHILE((uplimit-lowlimit)>bisectiontolinearsearchthreshold)
                      midpoint=(uplimit+lowlimit)/2
                      IF(matrix%ROW_INDICES(midpoint)>row_index) THEN
                        uplimit=midpoint
                      ELSE
                        lowlimit=midpoint
                      ENDIF
                    ENDDO
                    DO k=lowlimit,uplimit
                      IF(matrix%ROW_INDICES(k)==row_index) THEN
                        location=k
                        lowlimit=k+1
                        EXIT
                      ENDIF
                    ENDDO !k
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !i
                ENDDO !j
              CASE(matrix_row_column_storage_type)
                DO i=1,SIZE(row_indices,1)
                  row_index=row_indices(i)
                  DO j=1,SIZE(column_indices,1)
                    column_index=column_indices(j)
                    foundrow=.false.
                    location=1
                    DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                      IF(matrix%ROW_INDICES(location)==row_index) THEN
                        DO WHILE(.NOT.foundcolumn.AND.location<=matrix%SIZE)
                          IF(matrix%COLUMN_INDICES(location)==column_index.AND.matrix%ROW_INDICES(location)==row_index) THEN
                            foundcolumn=.true.
                          ELSE IF(matrix%ROW_INDICES(location)/=row_index) THEN
                            location=matrix%SIZE+1
                          ELSE
                            location=location+1
                          ENDIF
                        ENDDO
                      ELSE
                        location=location+1
                      ENDIF
                    ENDDO
                    IF(.NOT.(foundrow.AND.foundcolumn)) location=0
                    IF(location==0) THEN
                      local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                        & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                      CALL flagerror(local_error,err,error,*999)
                    ELSE
                      matrix%DATA_L(location)=matrix%DATA_L(location).OR.values(i,j)
                    ENDIF
                  ENDDO !j
                ENDDO !i
              CASE DEFAULT
                local_error="The matrix storage type of "//trim(number_to_vstring(matrix%STORAGE_TYPE,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_ADD_L2")
    RETURN
999 errorsexits("MATRIX_VALUES_ADD_L2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_add_l2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_intg(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_INTG",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  values(k)=0
                ELSE
                  values(k)=matrix%DATA_INTG(location)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_INTG")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_INTG",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_intg1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    INTEGER(INTG), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_INTG1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            VALUE=0
          ELSE
            VALUE=matrix%DATA_INTG(location)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_INTG1")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_INTG1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_intg1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_intg2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(OUT) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_INTG2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    values(i,j)=0
                  ELSE
                    values(i,j)=matrix%DATA_INTG(location)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_VALUES_GET_INTG2")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_INTG2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_intg2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_sp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_SP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  values(k)=0.0_sp
                ELSE
                  values(k)=matrix%DATA_SP(location)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_SP")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_SP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_sp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(SP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_SP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            VALUE=0.0_sp
          ELSE
            VALUE=matrix%DATA_SP(location)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_SP1")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_SP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_sp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_sp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(OUT) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_SP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    values(i,j)=0.0_sp
                  ELSE
                    values(i,j)=matrix%DATA_SP(location)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_SP2")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_SP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_sp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_dp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_DP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  values(k)=0.0_dp
                ELSE
                  values(k)=matrix%DATA_DP(location)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_DP")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_DP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_dp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_DP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            VALUE=0.0_dp
          ELSE
            VALUE=matrix%DATA_DP(location)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_DP1")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_DP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_dp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_dp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(OUT) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_DP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    values(i,j)=0.0_dp
                  ELSE
                    values(i,j)=matrix%DATA_DP(location)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_DP2")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_DP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_dp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_l(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_L",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  values(k)=.false.
                ELSE
                  values(k)=matrix%DATA_L(location)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_L")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_L",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_l

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_l1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    LOGICAL, INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_L1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            VALUE=.false.
          ELSE
            VALUE=matrix%DATA_L(location)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_L1")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_L1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_l1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_get_l2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(OUT) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_GET_L2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    values(i,j)=.false.
                  ELSE
                    values(i,j)=matrix%DATA_L(location)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_GET_L2")
    RETURN
999 errorsexits("MATRIX_VALUES_GET_L2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_get_l2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_intg(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_INTG",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_INTG(location)=values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_INTG")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_INTG",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_intg1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    INTEGER(INTG), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_INTG1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_INTG(location)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_INTG1")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_INTG1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_intg1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_intg2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_INTG2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_intg_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                      & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                    CALL flagerror(local_error,err,error,*999)
                  ELSE
                    matrix%DATA_INTG(location)=values(i,j)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the integer data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_VALUES_SET_INTG2")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_INTG2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_intg2
  
  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_sp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_SP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_SP(location)=values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_SP")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_SP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_sp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(SP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_SP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_SP(location)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_SP1")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_SP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_sp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_sp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(SP), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_SP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_sp_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                      & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                    CALL flagerror(local_error,err,error,*999)
                  ELSE
                    matrix%DATA_SP(location)=values(i,j)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the single precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF
    
    exits("MATRIX_VALUES_SET_SP2")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_SP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_sp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_dp(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_DP",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_DP(location)=values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_DP")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_DP",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_dp1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_DP1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_DP(location)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_DP1")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_DP1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_dp1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_dp2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    REAL(DP), INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_DP2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_dp_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(column_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                      & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                    CALL flagerror(local_error,err,error,*999)
                  ELSE
                    matrix%DATA_DP(location)=values(i,j)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the double precision data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_DP2")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_DP2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_dp2

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_l(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: k,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_L",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,1)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
              DO k=1,SIZE(row_indices,1)
                CALL matrix_storage_location_find(matrix,row_indices(k),column_indices(k),location,err,error,*999)
                IF(location==0) THEN
                  local_error="Row "//trim(number_to_vstring(row_indices(k),"*",err,error))//" and column "// &
                    & trim(number_to_vstring(column_indices(k),"*",err,error))//" does not exist in the matrix."
                  CALL flagerror(local_error,err,error,*999)
                ELSE
                  matrix%DATA_L(location)=values(k)
                ENDIF
              ENDDO !k
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_L")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_L",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_l

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_l1(MATRIX,ROW_INDEX,COLUMN_INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDEX
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDEX
    LOGICAL, INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_L1",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
          CALL matrix_storage_location_find(matrix,row_index,column_index,location,err,error,*999)
          IF(location==0) THEN
            local_error="Row "//trim(number_to_vstring(row_index,"*",err,error))//" and column "// &
              & trim(number_to_vstring(column_index,"*",err,error))//" does not exist in the matrix."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            matrix%DATA_L(location)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_L1")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_L1",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_l1

  !
  !================================================================================================================================
  !

  SUBROUTINE matrix_values_set_l2(MATRIX,ROW_INDICES,COLUMN_INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(matrix_type), POINTER :: MATRIX
    INTEGER(INTG), INTENT(IN) :: ROW_INDICES(:)
    INTEGER(INTG), INTENT(IN) :: COLUMN_INDICES(:)
    LOGICAL, INTENT(IN) :: VALUES(:,:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,j,LOCATION
    TYPE(varying_string) :: LOCAL_ERROR

    enters("MATRIX_VALUES_SET_L2",err,error,*999)

    IF(ASSOCIATED(matrix)) THEN
      IF(matrix%MATRIX_FINISHED) THEN
        IF(SIZE(row_indices,1)==SIZE(values,1)) THEN
          IF(SIZE(column_indices,1)==SIZE(values,2)) THEN
            IF(matrix%DATA_TYPE==matrix_vector_l_type) THEN
              DO i=1,SIZE(row_indices,1)
                DO j=1,SIZE(row_indices,1)
                  CALL matrix_storage_location_find(matrix,row_indices(i),column_indices(j),location,err,error,*999)
                  IF(location==0) THEN
                    local_error="Row "//trim(number_to_vstring(row_indices(i),"*",err,error))//" and column "// &
                      & trim(number_to_vstring(column_indices(j),"*",err,error))//" does not exist in the matrix."
                    CALL flagerror(local_error,err,error,*999)
                  ELSE
                    matrix%DATA_L(location)=values(i,j)
                  ENDIF
                ENDDO !j
              ENDDO !i
            ELSE
              local_error="The data type of "//trim(number_to_vstring(matrix%DATA_TYPE,"*",err,error))// &
                & " does not correspond to the logical data type of the given values."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ELSE
            local_error="The size of the column indices array ("// &
              & trim(number_to_vstring(SIZE(column_indices,1),"*",err,error))// &
              & ") does not conform to the number of columns in the values array ("// &
              & trim(number_to_vstring(SIZE(values,2),"*",err,error))//")."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the row indices array ("// &
            & trim(number_to_vstring(SIZE(row_indices,1),"*",err,error))// &
            & ") does not conform to the number of rows in the values array ("// &
            & trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The matrix has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Matrix is not associated.",err,error,*999)
    ENDIF

    exits("MATRIX_VALUES_SET_L2")
    RETURN
999 errorsexits("MATRIX_VALUES_SET_L2",err,error)
    RETURN 1
  END SUBROUTINE matrix_values_set_l2

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_all_values_set_intg(VECTOR,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_ALL_VALUES_SET_INTG",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
          vector%DATA_INTG=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_ALL_VALUES_SET_INTG")
    RETURN
999 errorsexits("VECTOR_ALL_VALUES_SET_INTG",err,error)
    RETURN 1
  END SUBROUTINE vector_all_values_set_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_all_values_set_sp(VECTOR,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    REAL(SP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_ALL_VALUES_SET_SP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
          vector%DATA_SP=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_ALL_VALUES_SET_SP")
    RETURN
999 errorsexits("VECTOR_ALL_VALUES_SET_SP",err,error)
    RETURN 1
  END SUBROUTINE vector_all_values_set_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_all_values_set_dp(VECTOR,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_ALL_VALUES_SET_DP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
          vector%DATA_DP=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_ALL_VALUES_SET_DP")
    RETURN
999 errorsexits("VECTOR_ALL_VALUES_SET_DP",err,error)
    RETURN 1
  END SUBROUTINE vector_all_values_set_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_all_values_set_l(VECTOR,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    LOGICAL, INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_ALL_VALUES_SET_L",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
          vector%DATA_L=VALUE
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_ALL_VALUES_SET_L")
    RETURN
999 errorsexits("VECTOR_ALL_VALUES_SET_L",err,error)
    RETURN 1
  END SUBROUTINE vector_all_values_set_l

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_create_finish(VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_CREATE_FINISH",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        CALL flagerror("Vector has been finished.",err,error,*999)
      ELSE
        IF(vector%SIZE>0) THEN
          SELECT CASE(vector%DATA_TYPE)
          CASE(matrix_vector_intg_type)
            ALLOCATE(vector%DATA_INTG(vector%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate vector integer data.",err,error,*999)
          CASE(matrix_vector_sp_type)
            ALLOCATE(vector%DATA_SP(vector%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate vector single precision data.",err,error,*999)
          CASE(matrix_vector_dp_type)
            ALLOCATE(vector%DATA_DP(vector%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate vector double precision data.",err,error,*999)
          CASE(matrix_vector_l_type)
            ALLOCATE(vector%DATA_L(vector%SIZE),stat=err)
            IF(err/=0) CALL flagerror("Could not allocate vector logical data.",err,error,*999)
          CASE DEFAULT
            local_error="The vector data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))//" is invalid."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        ENDIF
        vector%ID=matrix_vector_id
        matrix_vector_id=matrix_vector_id+1
        vector%VECTOR_FINISHED=.true.
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF
    
    exits("VECTOR_CREATE_FINISH")
    RETURN
999 errorsexits("VECTOR_CREATE_FINISH",err,error)
    RETURN 1
  END SUBROUTINE vector_create_finish

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_create_start(VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("VECTOR_CREATE_START",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      CALL flagerror("Vector is already associated.",err,error,*998)
    ELSE
      ALLOCATE(vector,stat=err)
      IF(err/=0) CALL flagerror("Could not allocate the vector.",err,error,*999)
      CALL vector_initialise(vector,err,error,*999)
      !Set the defaults
      vector%DATA_TYPE=matrix_vector_dp_type
    ENDIF
 
    exits("VECTOR_CREATE_START")
    RETURN
999 IF(ASSOCIATED(vector)) CALL vector_finalise(vector,err,error,*998)
998 errorsexits("VECTOR_CREATE_START",err,error)
    RETURN 1
  END SUBROUTINE vector_create_start

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_data_get_intg(VECTOR,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_DATA_GET_INTG",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(vector%VECTOR_FINISHED) THEN
          IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
            data=>vector%DATA_INTG
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the integer data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The vector has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF
    
    exits("VECTOR_DATA_GET_INTG")
    RETURN
999 errorsexits("VECTOR_DATA_GET_INTG",err,error)
    RETURN 1
  END SUBROUTINE vector_data_get_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_data_get_sp(VECTOR,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    REAL(SP), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_DATA_GET_SP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(vector%VECTOR_FINISHED) THEN
          IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
            data=>vector%DATA_SP
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the single precision data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The vector has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF
    
    exits("VECTOR_DATA_GET_SP")
    RETURN
999 errorsexits("VECTOR_DATA_GET_SP",err,error)
    RETURN 1
  END SUBROUTINE vector_data_get_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_data_get_dp(VECTOR,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    REAL(DP), POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_DATA_GET_DP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(vector%VECTOR_FINISHED) THEN
          IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
            data=>vector%DATA_DP
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the double precision data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The vector has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF
    
    exits("VECTOR_DATA_GET_DP")
    RETURN
999 errorsexits("VECTOR_DATA_GET_DP",err,error)
    RETURN 1
  END SUBROUTINE vector_data_get_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_data_get_l(VECTOR,DATA,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    LOGICAL, POINTER :: DATA(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_DATA_GET_L",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(ASSOCIATED(data)) THEN
        CALL flagerror("Data is already associated.",err,error,*999)
      ELSE
        NULLIFY(data)
        IF(vector%VECTOR_FINISHED) THEN
          IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
            data=>vector%DATA_L
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the logical data type of the requested values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("The vector has not been finished.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF
    
    exits("VECTOR_DATA_GET_L")
    RETURN
999 errorsexits("VECTOR_DATA_GET_L",err,error)
    RETURN 1
  END SUBROUTINE vector_data_get_l

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_datatypeget(vector,dataType,err,error,*)

    !Argument variables
    TYPE(vector_type), POINTER :: vector
    INTEGER(INTG), INTENT(OUT) :: dataType
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error

    enters("Vector_DataTypeGet",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(.NOT.vector%vector_finished) THEN
        CALL flag_error("The vector has not been finished.",err,error,*999)
      ELSE
        datatype=vector%data_type
      END IF
    ELSE
      CALL flag_error("Vector is not associated.",err,error,*999)
    END IF

    exits("Vector_DataTypeGet")
    RETURN
999 errorsexits("Vector_DataTypeGet",err,error)
    RETURN 1
  END SUBROUTINE vector_datatypeget

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_data_type_set(VECTOR,DATA_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: DATA_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_DATA_TYPE_SET",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        CALL flagerror("The vector has been finished.",err,error,*999)
      ELSE
        SELECT CASE(data_type)
        CASE(matrix_vector_intg_type)
          vector%DATA_TYPE=matrix_vector_intg_type
        CASE(matrix_vector_sp_type)
          vector%DATA_TYPE=matrix_vector_sp_type
        CASE(matrix_vector_dp_type)
          vector%DATA_TYPE=matrix_vector_dp_type
        CASE(matrix_vector_l_type)
          vector%DATA_TYPE=matrix_vector_l_type
        CASE DEFAULT
          local_error="The vector data type of "//trim(number_to_vstring(data_type,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_DATA_TYPE_SET")
    RETURN
999 errorsexits("VECTOR_DATA_TYPE_SET",err,error)
    RETURN 1
  END SUBROUTINE vector_data_type_set

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_destroy(VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("VECTOR_DESTROY",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      CALL vector_finalise(vector,err,error,*999)
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_DESTROY")
    RETURN
999 errorsexits("VECTOR_DESTROY",err,error)
    RETURN 1
  END SUBROUTINE vector_destroy

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_duplicate(VECTOR,NEW_VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    TYPE(vector_type), POINTER :: NEW_VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("VECTOR_DUPLICATE",err,error,*998)

    IF(ASSOCIATED(vector)) THEN
      IF(ASSOCIATED(new_vector)) THEN
        CALL flagerror("New vector is already associated.",err,error,*998)
      ELSE
        CALL vector_create_start(new_vector,err,error,*999)
        CALL vector_data_type_set(new_vector,vector%DATA_TYPE,err,error,*999)
        CALL vector_size_set(new_vector,vector%N,err,error,*999)
        CALL vector_create_finish(new_vector,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*998)
    ENDIF

    exits("VECTOR_DUPLICATE")
    RETURN
999 CALL vector_finalise(new_vector,err,error,*998)
998 errorsexits("VECTOR_DUPLICATE",err,error)
    RETURN 1
  END SUBROUTINE vector_duplicate

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_finalise(VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("VECTOR_FINALISE",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(ALLOCATED(vector%DATA_INTG)) DEALLOCATE(vector%DATA_INTG)
      IF(ALLOCATED(vector%DATA_SP)) DEALLOCATE(vector%DATA_SP)
      IF(ALLOCATED(vector%DATA_DP)) DEALLOCATE(vector%DATA_DP)
      IF(ALLOCATED(vector%DATA_L)) DEALLOCATE(vector%DATA_L)
      DEALLOCATE(vector)
    ENDIF

    exits("VECTOR_FINALISE")
    RETURN
999 errorsexits("VECTOR_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE vector_finalise

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_initialise(VECTOR,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    enters("VECTOR_INITIALISE",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      !!TODO: have a vector user number etc.
      vector%ID=0
      vector%VECTOR_FINISHED=.false.
      vector%N=0
      vector%DATA_TYPE=0
      vector%SIZE=0      
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_INITIALISE")
    RETURN
999 errorsexits("VECTOR_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE vector_initialise

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_size_set(VECTOR,N,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: N
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_SIZE_SET",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        CALL flagerror("The vector has been finished.",err,error,*999)
      ELSE
        IF(n>0) THEN
          vector%N=n
        ELSE
          local_error="The size of the vector ("//trim(number_to_vstring(n,"*",err,error))// &
            & ") is invalid. The number must be >0."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_SIZE_SET")
    RETURN
999 errorsexits("VECTOR_SIZE_SET",err,error)
    RETURN 1
  END SUBROUTINE vector_size_set

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_intg(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    INTEGER(INTG), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_INTG",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                values(i)=vector%DATA_INTG(k)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the integer data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_INTG")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_INTG",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_intg

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_intg1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    INTEGER(INTG), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_INTG1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            VALUE=vector%DATA_INTG(index)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_INTG1")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_INTG1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_intg1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_sp(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    REAL(SP), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_SP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                values(i)=vector%DATA_SP(k)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the single precision data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_SP")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_SP",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_sp1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    REAL(SP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_SP1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            VALUE=vector%DATA_SP(index)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_SP1")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_SP1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_sp1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_dp(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    REAL(DP), INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_DP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                values(i)=vector%DATA_DP(k)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the double precision data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_DP")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_DP",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_dp1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    REAL(DP), INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_DP1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            VALUE=vector%DATA_DP(index)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_DP1")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_DP1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_dp1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_l(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    LOGICAL, INTENT(OUT) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_L",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                values(i)=vector%DATA_L(k)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the logical data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")"
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_L")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_L",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_l

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_get_l1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    LOGICAL, INTENT(OUT) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_GET_L1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            VALUE=vector%DATA_L(index)
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_GET_L1")
    RETURN
999 errorsexits("VECTOR_VALUES_GET_L1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_get_l1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_intg(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    INTEGER(INTG), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_INTG",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                vector%DATA_INTG(k)=values(i)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the integer data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_INTG")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_INTG",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_intg

  !
  !================================================================================================================================
  !
  
  SUBROUTINE vector_values_set_intg1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    INTEGER(INTG), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_INTG1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_intg_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            vector%DATA_INTG(index)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the integer data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_INTG1")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_INTG1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_intg1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_sp(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    REAL(SP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_SP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
               vector%DATA_SP(k)=values(i)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the single precision data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_SP")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_SP",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_sp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_sp1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    REAL(SP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_SP1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_sp_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            vector%DATA_SP(index)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the single precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_SP1")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_SP1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_sp1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_dp(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    REAL(DP), INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_DP",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                vector%DATA_DP(k)=values(i)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the double precision data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_DP")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_DP",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_dp

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_dp1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    REAL(DP), INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_DP1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_dp_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            vector%DATA_DP(index)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the double precision data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_DP1")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_DP1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_dp1

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_l(VECTOR,INDICES,VALUES,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDICES(:)
    LOGICAL, INTENT(IN) :: VALUES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    INTEGER(INTG) :: i,k
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_L",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(SIZE(indices,1)==SIZE(values,1)) THEN
          IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
            DO i=1,SIZE(indices,1)
              k=indices(i)
              IF(k<1.OR.k>vector%N) THEN
                local_error="Index number "//trim(number_to_vstring(i,"*",err,error))//" is invalid. The index is "// &
                    & trim(number_to_vstring(k,"*",err,error))//" and it must be between 1 and "// &
                    & trim(number_to_vstring(vector%N,"*",err,error))//"."
                CALL flagerror(local_error,err,error,*999)
              ELSE
                vector%DATA_L(k)=values(i)
              ENDIF
            ENDDO !i
          ELSE
            local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
              & " does not correspond to the logical data type of the given values."
            CALL flagerror(local_error,err,error,*999)
          ENDIF
        ELSE
          local_error="The size of the indices array ("//trim(number_to_vstring(SIZE(indices,1),"*",err,error))// &
            & ") does not conform to the size of the values array ("//trim(number_to_vstring(SIZE(values,1),"*",err,error))//")."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_L")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_L",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_l

  !
  !================================================================================================================================
  !

  SUBROUTINE vector_values_set_l1(VECTOR,INDEX,VALUE,ERR,ERROR,*)

    !Argument variables
    TYPE(vector_type), POINTER :: VECTOR
    INTEGER(INTG), INTENT(IN) :: INDEX
    LOGICAL, INTENT(IN) :: VALUE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local variables
    TYPE(varying_string) :: LOCAL_ERROR

    enters("VECTOR_VALUES_SET_L1",err,error,*999)

    IF(ASSOCIATED(vector)) THEN
      IF(vector%VECTOR_FINISHED) THEN
        IF(vector%DATA_TYPE==matrix_vector_l_type) THEN
          IF(index<1.OR.index>vector%N) THEN
            local_error="The index value of "//trim(number_to_vstring(index,"*",err,error))// &
              & " is invalid. The index must be between 1 and  "//trim(number_to_vstring(vector%N,"*",err,error))//"."
            CALL flagerror(local_error,err,error,*999)
          ELSE
            vector%DATA_L(index)=VALUE
          ENDIF
        ELSE
          local_error="The data type of "//trim(number_to_vstring(vector%DATA_TYPE,"*",err,error))// &
            & " does not correspond to the logical data type of the given value."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("The vector has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Vector is not associated.",err,error,*999)
    ENDIF

    exits("VECTOR_VALUES_SET_L1")
    RETURN
999 errorsexits("VECTOR_VALUES_SET_L1",err,error)
    RETURN 1
  END SUBROUTINE vector_values_set_l1

  !
  !================================================================================================================================
  !
  
END MODULE matrix_vector