MatrixTypes_Native.f90 71.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief Futility-native implementations of MatrixTypes
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE MatrixTypes_Native
12
#include "Futility_DBC.h"
Graham, Aaron's avatar
Graham, Aaron committed
13
14
15
16
17
18
USE Futility_DBC
USE IntrType
USE ExceptionHandler
USE ParameterLists
USE Allocs
USE MatrixTypes_Base
19
20
USE ParallelEnv
USE Sorting
qicangsh's avatar
qicangsh committed
21

Graham, Aaron's avatar
Graham, Aaron committed
22
IMPLICIT NONE
23

24
25
26
27
#ifdef HAVE_MPI
#include <mpif.h>
#endif

Graham, Aaron's avatar
Graham, Aaron committed
28
PRIVATE
29
30
!
! List of public members
Graham, Aaron's avatar
Graham, Aaron committed
31
32
33
PUBLIC :: DenseSquareMatrixType
PUBLIC :: DenseRectMatrixType
PUBLIC :: TriDiagMatrixType
34
35
36
PUBLIC :: BandedMatrixType
PUBLIC :: DistributedBandedMatrixType
PUBLIC :: DistributedBlockBandedMatrixType
Graham, Aaron's avatar
Graham, Aaron committed
37
38
39
40
41
42
43
44
45
46
47
48
PUBLIC :: SparseMatrixType

!> @brief The extended type for dense square matrices
!>
!> This does not add a significant functionality over the intrinsic
!> allocatable arrays that are part of the Fortran language. It
!> is provided so as to be able to use the BLAS interfaces adapted
!> for the matrix type and it also gains some value by having the
!> isSymmetric and isInit attributes.
TYPE,EXTENDS(SquareMatrixType) :: DenseSquareMatrixType
  !> The values of the matrix
  REAL(SRK),ALLOCATABLE :: a(:,:)
49
50
!
!List of Type Bound Procedures
Graham, Aaron's avatar
Graham, Aaron committed
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
  CONTAINS
    !> @copybrief MatrixTypes::clear_DenseSquareMatrixType
    !> @copydetails MatrixTypes::clear_DenseSquareMatrixType
    PROCEDURE,PASS :: clear => clear_DenseSquareMatrixType
    !> @copybrief MatrixTypes::clear_DenseSquareMatrixType
    !> @copydetails MatrixTypes::clear_DenseSquareMatrixType
    PROCEDURE,PASS :: init => init_DenseSquareMatrixParam
    !> @copybrief MatrixTypes::set_DenseSquareMatrixType
    !> @copydetails MatrixTypes::set_DenseSquareMatrixType
    PROCEDURE,PASS :: set => set_DenseSquareMatrixType
    !> @copybrief MatrixTypes::get_DenseSquareMatrixType
    !> @copydetails MatrixTypes::get_DenseSquareMatrixType
    PROCEDURE,PASS :: get => get_DenseSquareMatrixType
    !> @copybrief MatrixTypes::transpose_DenseSquareMatrixType
    !> @copydetails MatrixTypes::transpose_DenseSquareMatrixType
    PROCEDURE,PASS :: transpose => transpose_DenseSquareMatrixType
    !> @copybrief MatrixTypes::zeroentries_DenseSquareMatrixType
    !> @copydetails MatrixTypes::zeroentries_DenseSquareMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_DenseSquareMatrixType
ENDTYPE DenseSquareMatrixType

!> @brief The extended type for dense rectangular matrices
TYPE,EXTENDS(RectMatrixType) :: DenseRectMatrixType
  !> The values of the matrix
  REAL(SRK),ALLOCATABLE :: a(:,:)
76
77
!
!List of Type Bound Procedures
Graham, Aaron's avatar
Graham, Aaron committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
  CONTAINS
    !> @copybrief MatrixTypes::clear_DenseRectMatrixType
    !> @copydetails MatrixTypes::clear_DenseRectMatrixType
    PROCEDURE,PASS :: clear => clear_DenseRectMatrixType
    !> @copybrief MatrixTypes::init_DenseRectMatrixType
    !> @copydetails MatrixTypes::init_DenseRectMatrixType
    PROCEDURE,PASS :: init => init_DenseRectMatrixParam
    !> @copybrief MatrixTypes::set_DenseRectMatrixType
    !> @copydetails MatrixTypes::set_DenseRectMatrixType
    PROCEDURE,PASS :: set => set_DenseRectMatrixType
    !> @copybrief MatrixTypes::get_DenseRectMatrixType
    !> @copydetails MatrixTypes::get_DenseRectMatrixType
    PROCEDURE,PASS :: get => get_DenseRectMatrixType
    !> @copybrief MatrixTypes::transpose_DenseRectMatrixType
    !> @copydetails MatrixTypes::transpose_DenseRectMatrixType
    PROCEDURE,PASS :: transpose => transpose_DenseRectMatrixType
    !> @copybrief MatrixTypes::zeroentries_DenseRectMatrixType
    !> @copydetails MatrixTypes::zeroentries_DenseRectMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_DenseRectMatrixType
ENDTYPE DenseRectMatrixType

!I think this may need to be revisited
!> @brief The extended type for tri-diagonal square matrices
TYPE,EXTENDS(SquareMatrixType) :: TriDiagMatrixType
  !> The values of the matrix
  REAL(SRK),ALLOCATABLE :: a(:,:)
104
105
!
!List of Type Bound Procedures
Graham, Aaron's avatar
Graham, Aaron committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
  CONTAINS
    !> @copybrief MatrixTypes::clear_TriDiagMatrixType
    !> @copydetails MatrixTypes::clear_TriDiagMatrixType
    PROCEDURE,PASS :: clear => clear_TriDiagMatrixType
    !> @copybrief MatrixTypes::init_TriDiagMatrixType
    !> @copydetails MatrixTypes::init_TriDiagMatrixType
    PROCEDURE,PASS :: init => init_TriDiagMatrixParam
    !> @copybrief MatrixTypes::set_TriDiagMatrixType
    !> @copydetails MatrixTypes::set_TriDiagMatrixType
    PROCEDURE,PASS :: set => set_TriDiagMatrixType
    !> @copybrief MatrixTypes::get_TriDiagMatrixType
    !> @copydetails MatrixTypes::get_TriDiagMatrixType
    PROCEDURE,PASS :: get => get_TriDiagMatrixType
    !> @copybrief MatrixTypes::transpose_TriDiagMatrixType
    !> @copydetails MatrixTypes::transpose_TriDiagMatrixType
    PROCEDURE,PASS :: transpose => transpose_TriDiagMatrixType
    !> @copybrief MatrixTypes::zeroentries_TriDiagMatrixType
    !> @copydetails MatrixTypes::zeroentries_TriDiagMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_TriDiagMatrixType
ENDTYPE TriDiagMatrixType

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
!> @brief Type used to hold the bands in the banded type
TYPE Band
  ! jIdx stores j index of each element in band
  INTEGER(SIK), ALLOCATABLE :: jIdx(:)
  ! elem stores the data in the band (compressed)
  REAL(SRK), ALLOCATABLE :: elem(:)
ENDTYPE Band

!> @brief The basic banded matrix type
TYPE,EXTENDS(MatrixType) :: BandedMatrixType
  !> Map of band indices stored (-m to n)
  INTEGER(SIK),ALLOCATABLE :: bandIdx(:)
  !> The bands stored in the matrix
  TYPE(Band),ALLOCATABLE :: bands(:)
  !> Number of nonzero elements
  INTEGER(SIK) :: nnz
  !> Number of columns
  INTEGER(SIK) :: m
  !> Counter to keep track of added elements before assembly
  INTEGER(SIK) :: counter
  !> Temporary containers used before (and deallocated after) assembly
  INTEGER(SIK), ALLOCATABLE :: iTmp(:),jTmp(:)
  REAL(SRK),ALLOCATABLE :: elemTmp(:)
  !> Logical flag for whether the matrix is assembled
  LOGICAL(SBK) :: isAssembled
  !> Logical flag for whether bands are stored in ascending or desending
  !> order by index
  LOGICAL(SBK) :: isReversed

!
!List of Type Bound Procedures
  CONTAINS
    !> @copybrief MatrixTypes::clear_BandedMatrixType
    !> @copydetails MatrixTypes::clear_BandedMatrixType
    PROCEDURE,PASS :: clear => clear_BandedMatrixType
    !> @copybrief MatrixTypes::init_BandedMatrixType
    !> @copydetails MatrixTypes::init_BandedMatrixType
    PROCEDURE,PASS :: init => init_BandedMatrixParam
    !> @copybrief MatrixTypes::assemble_BandedMatrixType
    !> @copydetails MatrixTypes::assemble_BandedMatrixType
    PROCEDURE,PASS :: assemble => assemble_BandedMatrixType
    !> @copybrief MatrixTypes::set_BandedMatrixType
    !> @copydetails MatrixTypes::set_BandedMatrixType
    PROCEDURE,PASS :: set => set_BandedMatrixType
    !> @copybrief MatrixTypes::get_BandedMatrixType
    !> @copydetails MatrixTypes::get_BandedMatrixType
    PROCEDURE,PASS :: get => get_BandedMatrixType
    !> @copybrief MatrixTypes::transpose_BandedMatrixType
    !> @copydetails MatrixTypes::transpose_BandedMatrixType
    PROCEDURE,PASS :: transpose => transpose_BandedMatrixType
    !> @copybrief MatrixTypes::zeroentries_BandedMatrixType
    !> @copydetails MatrixTypes::zeroentries_BandedMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_BandedMatrixType
    !> @copybrief MatrixTypes::binarySearch_BandedMatrixType
    !> @copydetails MatrixTypes::binarySearch_BandedMatrixType
    PROCEDURE,PASS,PRIVATE :: binarySearch => binarySearch_BandedMatrixType
ENDTYPE BandedMatrixType

!> @brief The distributed banded matrix type
TYPE,EXTENDS(DistributedMatrixType) :: DistributedBandedMatrixType
  !> Map of band indices stored (-m to n)
  INTEGER(SIK),ALLOCATABLE :: iOffsets(:),jOffsets(:)
  !> The column of banded matrix 'chunks' stored locally
  TYPE(BandedMatrixType),ALLOCATABLE :: chunks(:)
  !> Number of nonzero elements
192
  INTEGER(SLK) :: nnz
193
194
195
196
197
  !> Number of columns
  INTEGER(SIK) :: m
  !> Block size (smallest indivisble unit)
  INTEGER(SIK) :: blockSize
  !> Array of band sizes used to determine optimal communication
198
199
200
  INTEGER(SIK), ALLOCATABLE :: incIdxMap(:)
  INTEGER(SIK), ALLOCATABLE :: incIdxStt(:)
  INTEGER(SIK), ALLOCATABLE :: incIdxStp(:)
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  !> Temporary containers used before (and deallocated after) assembly
  INTEGER(SIK), ALLOCATABLE :: iTmp(:),jTmp(:)
  REAL(SRK),ALLOCATABLE :: elemTmp(:)
!
!List of Type Bound Procedures
  CONTAINS
    !> @copybrief MatrixTypes::clear_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::clear_DistributedBandedMatrixType
    PROCEDURE,PASS :: clear => clear_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::init_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::init_DistributedBandedMatrixType
    PROCEDURE,PASS :: init => init_DistributedBandedMatrixParam
    !> @copybrief MatrixTypes::assemble_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::assemble_DistributedBandedMatrixType
    PROCEDURE,PASS :: assemble => assemble_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::setrow_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::setrow_DistributedBandedMatrixType
    PROCEDURE,PASS :: setrow => setrow_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::set_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::set_DistributedBandedMatrixType
    PROCEDURE,PASS :: set => set_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::get_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::get_DistributedBandedMatrixType
    PROCEDURE,PASS :: get => get_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::transpose_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::transpose_DistributedBandedMatrixType
    PROCEDURE,PASS :: transpose => transpose_DistributedBandedMatrixType
    !> @copybrief MatrixTypes::zeroentries_DistributedBandedMatrixType
    !> @copydetails MatrixTypes::zeroentries_DistributedBandedMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_DistributedBandedMatrixType
ENDTYPE DistributedBandedMatrixType

!> @brief The block banded matrix type (designed for CMFD)
!> for now, the blocks will be assumed to be only on main diag
TYPE,EXTENDS(DistributedBandedMatrixType) :: DistributedBlockBandedMatrixType
  !> The block size is set by parent (DistrBanded::blockSize)
  !> The number of local blocks
  INTEGER(SIK) :: nlocalBlocks
  !> Block offset for this processor
  INTEGER(SIK) :: blockOffset
  !> Mask to effectively zero out block vals
  LOGICAL(SBK) :: blockMask
  !> The dense block container
  TYPE(DenseSquareMatrixType),ALLOCATABLE :: blocks(:)
!
!List of Type Bound Procedures
  CONTAINS
    !> @copybrief MatrixTypes::clear_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::clear_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: clear => clear_DistributedBlockBandedMatrixType
    !> @copybrief MatrixTypes::init_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::init_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: init => init_DistributedBlockBandedMatrixParam
    ! The assemble routine will remain unchanged
    ! The setRow routine will remain unchanged (unimplemented)
    !> @copybrief MatrixTypes::set_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::set_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: set => set_DistributedBlockBandedMatrixType
    !> @copybrief MatrixTypes::get_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::get_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: get => get_DistributedBlockBandedMatrixType
    ! the transpose routine will remain unchanged (unimplemented)
    !> @copybrief MatrixTypes::zeroentries_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::zeroentries_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_DistributedBlockBandedMatrixType
    !> @copybrief MatrixTypes::setBlockMask_DistributedBlockBandedMatrixType
    !> @copydetails MatrixTypes::setBlockMask_DistributedBlockBandedMatrixType
    PROCEDURE,PASS :: setBlockMask => setBlockMask_DistributedBlockBandedMatrixType
ENDTYPE DistributedBlockBandedMatrixType


Graham, Aaron's avatar
Graham, Aaron committed
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
!> @brief The basic sparse matrix type
!>
!> Matrix uses compressed sparse row storage format,
!> as defined by Intel's MKL
TYPE,EXTENDS(MatrixType) :: SparseMatrixType
  !> Number of non-zero entries in the matrix
  INTEGER(SIK) :: nnz=0
  !> The number of elements in each row
  INTEGER(SIK),ALLOCATABLE :: ia(:)
  !> The column indices for each element
  INTEGER(SIK),ALLOCATABLE :: ja(:) !columns
  !> The values of the matrix
  REAL(SRK),ALLOCATABLE :: a(:) !values
  !> A counter for the current location in a(:) and ja(:)
  INTEGER(SIK) :: jCount=0
  !> A variable to store the previous row entered in set_shape
  INTEGER(SIK) :: iPrev=0
  !> A variable to store the previous column entered in set_shape
  INTEGER(SIK) :: jPrev=0
291
292
!
!List of Type Bound Procedures
Graham, Aaron's avatar
Graham, Aaron committed
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
  CONTAINS
    !> @copybrief MatrixTypes::clear_SparseMatrixType
    !> @copydetails MatrixTypes::clear_SparseMatrixType
    PROCEDURE,PASS :: clear => clear_SparseMatrixType
    !> @copybrief MatrixTypes::init_SparseMatrixType
    !> @copydetails MatrixTypes::init_SparseMatrixType
    PROCEDURE,PASS :: init => init_SparseMatrixParam
    !> @copybrief MatrixTypes::set_SparseMatrixType
    !> @copydetails MatrixTypes::set_SparseMatrixType
    PROCEDURE,PASS :: set => set_SparseMatrixType
    !> @copybrief MatrixTypes::setRow_SparseMatrixType
    !> @copydetails MatrixTypes::setRow_SparseMatrixType
    PROCEDURE,PASS :: setRow => setRow_SparseMatrixType
    !> @copybrief MatrixTypes::set_shape_SparseMatrixType
    !> @copydetails MatrixTypes::set_shape_SparseMatrixType
    PROCEDURE,PASS :: setShape => set_shape_SparseMatrixType
    !> @copybrief MatrixTypes::get_SparseMatrixType
    !> @copydetails MatrixTypes::get_SparseMatrixType
    PROCEDURE,PASS :: get => get_SparseMatrixType
    !> @copybrief MatrixTypes::transpose_SparseMatrixType
    !> @copydetails MatrixTypes::transpose_SparseMatrixType
    PROCEDURE,PASS :: transpose => transpose_SparseMatrixType
    !> @copybrief MatrixTypes::zeroentries_SparseMatrixType
    !> @copydetails MatrixTypes::zeroentries_SparseMatrixType
    PROCEDURE,PASS :: zeroentries => zeroentries_SparseMatrixType
ENDTYPE SparseMatrixType

!> Name of module
CHARACTER(LEN=*),PARAMETER :: modName='MATRIXTYPES_NATIVE'
322
323
324

!
!===============================================================================
Graham, Aaron's avatar
Graham, Aaron committed
325
CONTAINS
326
327
328
329
!
!-------------------------------------------------------------------------------
!> @brief Initializes Sparse Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
330
!> @param Params the parameter list
331
!>
Graham, Aaron's avatar
Graham, Aaron committed
332
333
334
335
336
SUBROUTINE init_SparseMatrixParam(matrix,Params)
  CLASS(SparseMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
  TYPE(ParamType) :: validParams
  INTEGER(SIK) :: n,nnz
337

338
339
  REQUIRE(.NOT.matrix%isInit)

Graham, Aaron's avatar
Graham, Aaron committed
340
341
  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()
342

Graham, Aaron's avatar
Graham, Aaron committed
343
344
345
  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(SparseMatrixType_reqParams)
346

Graham, Aaron's avatar
Graham, Aaron committed
347
348
349
350
  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
  CALL validParams%get('MatrixType->nnz',nnz)
  CALL validParams%clear()
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
  REQUIRE(n > 0)
  REQUIRE(nnz > 0)

  matrix%isInit=.TRUE.
  matrix%n=n
  matrix%nnz=nnz
  matrix%jCount=0
  matrix%iPrev=0
  matrix%jPrev=0
  !regardless of sparsity, SIZE(ia)=n+1
  CALL dmallocA(matrix%ia,matrix%n+1)
  CALL dmallocA(matrix%a,matrix%nnz)
  CALL dmallocA(matrix%ja,matrix%nnz)
  !last entry of ia is known in advanced
  !this is per the intel MKL format
  matrix%ia(matrix%n+1)=matrix%nnz+1
367

Graham, Aaron's avatar
Graham, Aaron committed
368
ENDSUBROUTINE init_SparseMatrixParam
369
370
371
372
!
!-------------------------------------------------------------------------------
!> @brief Initializes Tridiagonal Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
373
!> @param Params the parameter list
374
!>
Graham, Aaron's avatar
Graham, Aaron committed
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
SUBROUTINE init_TriDiagMatrixParam(matrix,Params)
  CHARACTER(LEN=*),PARAMETER :: myName='init_TriDiagMatrixParam'
  CLASS(TriDiagMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
  TYPE(ParamType) :: validParams
  INTEGER(SIK) :: n
  LOGICAL(SBK) :: isSym

  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()

  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(TriDiagMatrixType_reqParams)

  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
  CALL validParams%get('MatrixType->isSym',isSym)
  CALL validParams%clear()

  IF(.NOT. matrix%isInit) THEN
    IF(n < 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '// &
          modName//'::'//myName//' - Number of rows (n) must be '// &
          'greater than 1!')
    ELSE
      matrix%isInit=.TRUE.
      matrix%n=n
      IF(isSym) THEN
        matrix%isSymmetric=.TRUE.
405
      ELSE
Graham, Aaron's avatar
Graham, Aaron committed
406
        matrix%isSymmetric=.FALSE.
407
      ENDIF
Graham, Aaron's avatar
Graham, Aaron committed
408
409
410
411
412
413
414
      CALL dmallocA(matrix%a,3,n)
    ENDIF
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF
ENDSUBROUTINE init_TriDiagMatrixParam
415
416
!
!-------------------------------------------------------------------------------
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
!> @brief Initializes Banded Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
!> @param Params the parameter list
!>
SUBROUTINE init_BandedMatrixParam(matrix,Params)
  CHARACTER(LEN=*),PARAMETER :: myName='init_BandedMatrixParam'
  CLASS(BandedMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
  TYPE(ParamType) :: validParams
  INTEGER(SIK) :: n,m,nnz

  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()

  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(BandedMatrixType_reqParams)

  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
437
438
  m=n
  IF (validParams%has('MatrixType->m')) CALL validParams%get('MatrixType->m',m)
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
  CALL validParams%get('MatrixType->nnz',nnz)
  CALL validParams%clear()

  IF(.NOT. matrix%isInit) THEN
    IF(n <= 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of rows (n) must be greater than 1!')
    ELSEIF(m <= 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of columns (m) must be greater than 1!')
    ELSEIF(nnz < 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of nonzero elements (nnz) must be greater than 0!')
    ELSE
      ALLOCATE(matrix%iTmp(nnz))
      ALLOCATE(matrix%jTmp(nnz))
      ALLOCATE(matrix%elemTmp(nnz))

      matrix%isInit=.TRUE.
      matrix%isAssembled=.FaLSE.
      matrix%isReversed=.FALSE.
      matrix%counter=0_SRK
      matrix%n=n
      matrix%m=m
      matrix%nnz=nnz
    ENDIF
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF
ENDSUBROUTINE init_BandedMatrixParam
!
!-------------------------------------------------------------------------------
!> @brief Initializes Distributed Banded Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
!> @param Params the parameter list
!>
SUBROUTINE init_DistributedBandedMatrixParam(matrix,Params)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
#ifdef HAVE_MPI
  CHARACTER(LEN=*),PARAMETER :: myName='init_DistributedBandedMatrixParam'
  TYPE(ParamType) :: validParams
482
483
484
  INTEGER(SIK),ALLOCATABLE :: onnz(:),dnnz(:)
  INTEGER(SIK) :: n,m,commID,rank,mpierr,nproc,i,blocksize,nlocal,nnz_int
  INTEGER(SLK) :: nnz
485
486
487
488
489
490
491
492
493
494
495
496
497

  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()

  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(DistributedBandedMatrixType_reqParams, &
      DistributedBandedMatrixType_optParams)

  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
  CALL validParams%get('MatrixType->m',m)
  CALL validParams%get('MatrixType->MPI_Comm_ID',commID)
498
  CALL validParams%get('MatrixType->nnz',nnz_int)
499
500
  CALL validParams%get('MatrixType->blockSize',blockSize)
  CALL validParams%get('MatrixType->nlocal',nlocal)
501
502
503
504
505
506
507
508
509
510
511
512
513
  nnz=INT(nnz_int,KIND=8)
  IF(nnz_int <= 1) THEN
    IF (validParams%has("MatrixType->dnnz") .AND. validParams%has("MatrixType->onnz")) THEN
      nnz=0_SLK
      CALL validParams%get("MatrixType->onnz",onnz)
      CALL validParams%get("MatrixType->dnnz",dnnz)
      DO i=1,SIZE(onnz)
        nnz=nnz+INT(onnz(i),KIND=SLK)
        nnz=nnz+INT(dnnz(i),KIND=SLK)
      ENDDO
      CALL MPI_AllReduce(nnz,MPI_IN_PLACE,1,MPI_LONG,MPI_SUM,commID,mpierr)
    ENDIF
  ENDIF
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548

  IF(.NOT. matrix%isInit) THEN
    IF(n <= 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of rows (n) must be greater than 1!')
    ELSEIF(m <= 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of columns (m) must be greater than 1!')
    ELSEIF(nnz <= 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - Number of nonzero entries (nnz) must be greater than 0!')
    ELSEIF(commID == MPI_COMM_NULL) THEN
      CALL eMatrixType%raiseError('Incorrect input to '//modName//'::'//myName// &
          ' - MPI communicator cannot have the same value as MPI_COMM_NULL')
    ELSE
      CALL MPI_Comm_rank(commID,rank,mpierr)
      CALL MPI_Comm_size(commID,nproc,mpierr)

      REQUIRE(MOD(n,blockSize)==0)
      REQUIRE(MOD(m,blockSize)==0)

      ALLOCATE(matrix%jOffsets(nproc+1))
      ALLOCATE(matrix%iOffsets(nproc+1))
      matrix%jOffsets(1) = 0
      matrix%iOffsets(1) = 0

      matrix%nlocal = 0
      n=n/blockSize
      m=m/blockSize
      ! We assume that the matrix is somewhat evenly balanced
      ! among its columns; that is, no one chunk of columns contains
      ! more than twice the average nonzero entries
      ALLOCATE(matrix%iTmp(2*nnz/nProc))
      ALLOCATE(matrix%jTmp(2*nnz/nProc))
      ALLOCATE(matrix%elemTmp(2*nnz/nProc))
549
      IF(nlocal < 0) THEN
550
551
552
        ! Automatic division of matrix by greedy partitioning of blocks
        DO i=2,nproc+1
          matrix%jOffsets(i)=matrix%jOffsets(i-1)+m/nproc
553
          IF(MOD(m,nproc) > i-2) THEN
554
555
556
557
558
559
            matrix%jOffsets(i)=matrix%jOffsets(i)+1
          ENDIF
        ENDDO

        DO i=2,nproc+1
          matrix%iOffsets(i)=matrix%iOffsets(i-1)+n/nproc
560
          IF(MOD(n,nproc) > i-2) THEN
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
            matrix%iOffsets(i)=matrix%iOffsets(i)+1
          ENDIF
        ENDDO
        ! Separate i/jOffsets are maintained in the case of a rectangular
        ! matrix
        matrix%iOffsets=matrix%iOffsets*blockSize
        matrix%jOffsets=matrix%jOffsets*blockSize
      ELSE
        ! Partitioning explicitly specified (square only)
        REQUIRE(MOD(nlocal,blockSize) == 0)
        REQUIRE(m == n)
        ! Gather all nlocal and sum going forward
        CALL MPI_AllGather(nlocal,1,MPI_INTEGER,matrix%iOffsets(2),1, &
            MPI_INTEGER,commID,mpierr)
        DO i=2,nproc+1
          matrix%iOffsets(i)=matrix%iOffsets(i)+matrix%iOffsets(i-1)
        ENDDO
        REQUIRE(matrix%iOffsets(SIZE(matrix%iOffsets)) == n*blockSize)
        matrix%jOffsets(:)=matrix%iOffsets(:)
      ENDIF

      matrix%isInit=.TRUE.
      matrix%isAssembled=.FALSE.
      matrix%n=n*blockSize
      matrix%m=m*blockSize
      matrix%nnz=nnz
      matrix%blockSize = blockSize
      matrix%comm=commID
    ENDIF
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF
594
  CALL validParams%clear()
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
#endif
ENDSUBROUTINE init_DistributedBandedMatrixParam

!
!-------------------------------------------------------------------------------
!> @brief Initializes Distributed Block-Banded Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
!> @param Params the parameter list
!>
SUBROUTINE init_DistributedBlockBandedMatrixParam(matrix,Params)
  CLASS(DistributedBlockBandedMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
#ifdef HAVE_MPI
  CHARACTER(LEN=*),PARAMETER :: myName='init_DistributedBlockBandedMatrixParam'
  TYPE(ParamType) :: validParams,blockParams
  INTEGER(SIK) :: n,rank,mpierr,commID,nproc,i,blockSize,nlocal
  INTEGER(SIK),ALLOCATABLE :: blkOffsetTmp(:)

613
  IF(.NOT. matrix%isInit) THEN
614
    !Check to set up required and optional param lists.
615
    IF(.NOT. MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643

    !Validate against the reqParams and OptParams
    validParams=Params
    CALL validParams%validate(DistributedBlockBandedMatrixType_reqParams, &
        DistributedBlockBandedMatrixType_optParams)

    ! Pull Data From Parameter List
    CALL validParams%get('MatrixType->n',n)
    CALL validParams%get('MatrixType->MPI_Comm_ID',commID)
    CALL validParams%get('MatrixType->blockSize',blockSize)
    CALL validParams%get('MatrixType->nlocal',nlocal)

    REQUIRE(.NOT. matrix%isInit)
    REQUIRE(blockSize > 0)
    REQUIRE(MOD(n,blockSize) == 0)

    matrix%blockMask=.FALSE.

    ! Allocate the blocks, then adjust the parameter list and proceed to
    ! %parent%

    ! Default to greedy partitioning, respecting blockSize
    CALL MPI_Comm_rank(commID,rank,mpierr)
    CALL MPI_Comm_size(commID,nproc,mpierr)
    REQUIRE(n/blockSize >= nproc)

    n=n/blockSize
    ! If nlocal < 0, default to greedy partitioning
644
645
    IF(nlocal < 0) THEN
      IF(rank < MOD(n,nproc)) THEN
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
        matrix%nlocalBlocks=n/nproc+1
        matrix%blockOffset=(rank)*(n/nproc+1)
      ELSE
        matrix%nlocalBlocks=n/nproc
        matrix%blockOffset=(rank)*(n/nproc)+MOD(n,nproc)
      ENDIF
    ELSE
      ALLOCATE(blkOffsetTmp(nproc))
      REQUIRE(MOD(nlocal,blockSize) == 0)
      matrix%nlocalBlocks=nlocal/blockSize
      CALL MPI_AllGather(matrix%nlocalBlocks,1,MPI_INTEGER,blkOffsetTmp,1, &
          MPI_INTEGER,commID,mpierr)
      matrix%blockOffset=SUM(blkOffsetTmp(1:rank))
      DEALLOCATE(blkOffsetTmp)
    ENDIF

    ALLOCATE(matrix%blocks(matrix%nlocalBlocks))
    CALL blockParams%clear()
    CALL blockParams%add('MatrixType->n',blockSize)
    CALL blockParams%add('MatrixType->isSym',.FALSE.)

    DO i=1,matrix%nlocalblocks
      CALL matrix%blocks(i)%init(blockParams)
    ENDDO
670
    CALL blockParams%clear()
671
672
673

    CALL validParams%add('MatrixType->m',n*blockSize)
    CALL matrix%DistributedBandedMatrixType%init(validParams)
674
    CALL validParams%clear()
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF

#endif
ENDSUBROUTINE init_DistributedBlockBandedMatrixParam

!
!-------------------------------------------------------------------------------
!> @brief Assembles Serial Banded Matrix Type. This rearranges values set into
!>        matrix bands, and sets the structure to read-only
!> @param matrix the matrix type to act on
!>
SUBROUTINE assemble_BandedMatrixType(thisMatrix)
  CLASS(BandedMatrixType),INTENT(INOUT) :: thisMatrix
  INTEGER(SIK),ALLOCATABLE :: numOnDiag(:),idxOrig(:)
  INTEGER(SIK),ALLOCATABLE :: iSort(:),jSort(:)
  REAL(SRK),ALLOCATABLE :: valSort(:)
  INTEGER(SIK) :: counter,nBands,i,j
  INTEGER(SLK) :: iLong,jLong,nLong,mLong
  INTEGER(SLK),ALLOCATABLE :: diagRank(:)

  REQUIRE(thisMatrix%isInit)

700
701
702
703
  ALLOCATE(diagRank(thisMatrix%counter))
  ALLOCATE(idxOrig(thisMatrix%counter))
  IF (thisMatrix%isAssembled) RETURN
  DO i=1,thisMatrix%counter
704
705
706
707
    iLong=INT(thisMatrix%iTmp(i),KIND=SLK)
    jLong=INT(thisMatrix%jTmp(i),KIND=SLK)
    nLong=INT(thisMatrix%n,KIND=SLK)
    mLong=INT(thisMatrix%m,KIND=SLK)
708
709
710
711
    idxOrig(i)=i
    ! The diagonal rank counts upwards as one moves southeast in the matrix,
    ! starting at a large negative number in the bottom left, and reaching
    ! zero near the middle
712
    IF(thisMatrix%n-thisMatrix%iTmp(i)+1+thisMatrix%jTmp(i) &
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
        <= MAX(thisMatrix%m,thisMatrix%n)) THEN
      diagRank(i)=-mLong*nLong/2+jLong-1+((nLong-iLong+jLong-1)*(nLong-iLong+jLong))/2
    ELSE
      diagRank(i) = nLong*mLong/2-(mLong-jLong+((iLong+mLong-jLong-1)*(iLong+mLong-jLong))/2)
    ENDIF
  ENDDO

  CALL sort(diagRank,idxOrig)
  ! Now update iTmp,jTmp,etc
  ALLOCATE(iSort(SIZE(idxOrig)))
  iSort = thisMatrix%iTmp(idxOrig)
  DEALLOCATE(thisMatrix%iTmp)
  jSort = thisMatrix%jTmp(idxOrig)
  DEALLOCATE(thisMatrix%jTmp)
  valSort = thisMatrix%elemTmp(idxOrig)
  DEALLOCATE(thisMatrix%elemTmp)

  ! Find number of elements in each band
  ALLOCATE(numOnDiag(thisMatrix%m+thisMatrix%n-1))
  numOnDiag=0_SIK
  DO i=1,SIZE(iSort)
    numOnDiag(jSort(i)-iSort(i)+thisMatrix%n)=1+numOnDiag(jSort(i)-iSort(i)+thisMatrix%n)
  ENDDO
  ! Now count up the number of bands necessary
  nBands=0
  DO i=1,SIZE(numOnDiag)
739
    IF(numOnDiag(i)/=0) THEN
740
741
742
743
744
745
746
747
      nBands=nBands+1
    ENDIF
  ENDDO
  ALLOCATE(thisMatrix%bandIdx(nBands))
  ALLOCATE(thisMatrix%bands(nBands))
  counter=1
  ! Allocate the bands
  DO i=1,SIZE(numOnDiag)
748
    IF(numOnDiag(i)/=0) THEN
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
      ALLOCATE(thisMatrix%bands(counter)%jIdx(numOnDiag(i)))
      ALLOCATE(thisMatrix%bands(counter)%elem(numOnDiag(i)))
      thisMatrix%bandIdx(counter)=i-thisMatrix%n
      counter=counter+1
    ENDIF
  ENDDO
  counter=1
  ! Set the bands
  DO i=1,SIZE(thisMatrix%bandIdx)
    DO j=1,SIZE(thisMatrix%bands(i)%jIdx)
      thisMatrix%bands(i)%jIdx(j)=jSort(counter)
      thisMatrix%bands(i)%elem(j)=valSort(counter)
      counter=counter+1
    ENDDO
  ENDDO
  thisMatrix%isAssembled=.TRUE.
ENDSUBROUTINE assemble_BandedMatrixType
!
!-------------------------------------------------------------------------------
!> @brief Assembles Distributed Banded Matrix Type. This sets some values and
!>        assembles all chunks, setting the structure to read-only
!> @param thisMatrix the matrix type to act on
!> @param ierr optional error code return (not used for native types)
!>
SUBROUTINE assemble_DistributedBandedMatrixType(thisMatrix,ierr)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: thisMatrix
  INTEGER(SIK),INTENT(OUT),OPTIONAL :: ierr
#ifdef HAVE_MPI
  TYPE(ParamType) :: bandedPList
778
779
  INTEGER(SIK) :: mpierr,rank,nproc,i,j,nTransmit,stt,stp
  INTEGER(SIK),ALLOCATABLE :: nnzPerChunk(:),nRecv(:),idxMapTmp(:)
780
781
782
783
784

  REQUIRE(thisMatrix%isInit)
  CALL MPI_Comm_rank(thisMatrix%comm,rank,mpierr)
  CALL MPI_Comm_size(thisMatrix%comm,nproc,mpierr)
  ALLOCATE(nnzPerChunk(nProc))
785
  ALLOCATE(nRecv(nProc))
786
787
788
  nnzPerChunk=0
  DO i=1,thisMatrix%nLocal
    DO j=1,nproc
789
      IF(thisMatrix%iTmp(i) > thisMatrix%iOffsets(j) .AND. &
790
791
792
793
794
795
796
797
798
799
          thisMatrix%iTmp(i) <= thisMatrix%iOffsets(j+1)) THEN
        nnzPerChunk(j)=nnzPerChunk(j)+1
        EXIT
      ENDIF
    ENDDO
  ENDDO
  CALL bandedPList%clear()

  ALLOCATE(thisMatrix%chunks(nProc))
  DO i=1,nProc
800
    IF(nnzPerChunk(i) > 0) THEN
801
802
803
804
805
806
807
808
      CALL bandedPList%add('MatrixType->nnz',nnzPerChunk(i))
      CALL bandedPList%add('MatrixType->n',thisMatrix%iOffsets(i+1) &
          -thisMatrix%iOffsets(i))
      CALL bandedPList%add('MatrixType->m',thisMatrix%jOffsets(rank+2) &
          -thisMatrix%jOffsets(rank+1))
      CALL thisMatrix%chunks(i)%init(bandedPList)

      DO j=1,thisMatrix%nLocal
809
        IF(thisMatrix%iTmp(j) > thisMatrix%iOffsets(i) .AND. &
810
811
812
813
814
815
816
817
            thisMatrix%iTmp(j) <= thisMatrix%iOffsets(i+1)) THEN
          CALL thisMatrix%chunks(i)%set(thisMatrix%iTmp(j)-thisMatrix%iOffsets(i), &
              thisMatrix%jTmp(j)-thisMatrix%jOffsets(rank+1),thisMatrix%elemTmp(j))
        ENDIF
      ENDDO
      CALL thisMatrix%chunks(i)%assemble()
      CALL bandedPList%clear()
    ENDIF
818
819
820
821
822
823
824
825
826

    ! Decide the sending method to use. 'Full' chunks will send their whole
    ! contents while 'sparse' chunks will send band products.
    ! Set nTransmit to the number of REAL's to be sent from that processor to
    ! this one.
    !  >0: This many
    !   0: None (chunk empty)
    !  -1: All
    nTransmit=0
827
    IF(thisMatrix%chunks(i)%isInit) THEN
828
829
830
831
      IF(thisMatrix%chunks(i)%nnz < thisMatrix%chunks(i)%n) THEN
        DO j=1,SIZE(thisMatrix%chunks(i)%bandIdx)
          nTransmit=nTransmit+SIZE(thisMatrix%chunks(i)%bands(j)%jIdx)
        ENDDO
832
      ELSE
833
        nTransmit=-1
834
835
      ENDIF
    ENDIF
836
    ! This info must be disseminated through all procs
837

838
839
    ! Gather to rank i-1 into nRecv(:)
    CALL MPI_Gather(nTransmit,1,MPI_INTEGER,nRecv,1,MPI_INTEGER,i-1, &
840
841
        thisMatrix%comm,mpierr)
    ! Have rank i-1 allocate
842
    IF(rank == i-1) THEN
843
844
845
846
847
848
849
850
851
852
853
854
855
856
      nRecv(i)=0
      nTransmit=0
      DO j=1,SIZE(nRecv)
        IF (nRecv(j) > 0) THEN
          nTransmit=nTransmit+nRecv(j)
        ENDIF
      ENDDO
      ALLOCATE(thisMatrix%incIdxMap(nTransmit))
      ALLOCATE(thisMatrix%incIdxStt(nProc))
      ALLOCATE(thisMatrix%incIdxStp(nProc))
      thisMatrix%incIdxStt=0
      thisMatrix%incIdxStp=0
      ! Use nTransmit as counter
      nTransmit=1
857
      DO j=1,nProc
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
        IF(j == i) THEN
          thisMatrix%incIdxStt(j)=0
          thisMatrix%incIdxStp(j)=0
        ELSE
          IF(nRecv(j) > 0) THEN
            thisMatrix%incIdxStt(j)=nTransmit
            thisMatrix%incIdxStp(j)=nTransmit+nRecv(j)-1
            ! We receive data (and don't want to bother sending idx's at mvmult
            ! time)
            CALL MPI_Recv(thisMatrix%incIdxMap(nTransmit:nTransmit+nRecv(j)-1), &
                nRecv(j),MPI_INTEGER,j-1,0,thisMatrix%comm,MPI_STATUS_IGNORE, mpierr)
            nTransmit=nTransmit+nRecv(j)
          ELSE IF (nRecv(j) == -1) THEN
            thisMatrix%incIdxStt(j)=-1
            thisMatrix%incIdxStp(j)=-1
873
874
875
876
          ENDIF
        ENDIF
      ENDDO
    ELSE
877
878
879
880
881
882
883
884
      ! Build the index map
      IF(nTransmit > 0) THEN
        ALLOCATE(idxMapTmp(nTransmit))
        stt=1
        DO j=1,SIZE(thisMatrix%chunks(i)%bandIdx)
          stp=stt+SIZE(thisMatrix%chunks(i)%bands(j)%jIdx)-1
          idxMapTmp(stt:stp)=thisMatrix%chunks(i)%bands(j)%jIdx-thisMatrix%chunks(i)%bandIdx(j)
          stt=stp+1
885
        ENDDO
886
        CALL MPI_Send(idxMapTmp,nTransmit,MPI_INTEGER,i-1,0, &
887
            thisMatrix%comm,mpierr)
888
        DEALLOCATE(idxMapTmp)
889
890
891
892
893
894
895
896
897
898
899
900
      ENDIF
    ENDIF
  ENDDO

  thisMatrix%isAssembled = .TRUE.
  DEALLOCATE(thisMatrix%iTmp)
  DEALLOCATE(thisMatrix%jTmp)
  DEALLOCATE(thisMatrix%elemTmp)
#endif
ENDSUBROUTINE assemble_DistributedBandedMatrixType
!
!-------------------------------------------------------------------------------
901
902
!> @brief Initializes Dense Rectangular Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
903
!> @param Params the parameter list
904
!>
Graham, Aaron's avatar
Graham, Aaron committed
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
SUBROUTINE init_DenseRectMatrixParam(matrix,Params)
  CHARACTER(LEN=*),PARAMETER :: myName='init_DenseRectMatrixParam'
  CLASS(DenseRectMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
  TYPE(ParamType) :: validParams
  INTEGER(SIK) :: n,m

  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()

  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(DenseRectMatrixType_reqParams)

  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
  CALL validParams%get('MatrixType->m',m)
  CALL validParams%clear()

  IF(.NOT. matrix%isInit) THEN
    IF(n < 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '// &
          modName//'::'//myName//' - Number of rows (n) must'// &
          ' be greater than 1!')
    ELSEIF(m < 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '// &
          modName//'::'//myName//' - Number of columns (m) must'// &
          ' be greater than 1!')
    ELSE
      matrix%isInit=.TRUE.
      matrix%n=n
      matrix%m=m
      CALL dmallocA(matrix%a,n,m)
    ENDIF
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF
ENDSUBROUTINE init_DenseRectMatrixParam
944
945
946
947
!
!-------------------------------------------------------------------------------
!> @brief Initializes Dense Square Matrix Type with a Parameter List
!> @param matrix the matrix type to act on
948
!> @param Params the parameter list
949
!>
Graham, Aaron's avatar
Graham, Aaron committed
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
SUBROUTINE init_DenseSquareMatrixParam(matrix,Params)
  CHARACTER(LEN=*),PARAMETER :: myName='init_DenseSquareMatrixParam'
  CLASS(DenseSquareMatrixType),INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: Params
  TYPE(ParamType) :: validParams
  INTEGER(SIK) :: n
  LOGICAL(SBK) :: isSym

  !Check to set up required and optional param lists.
  IF(.NOT.MatrixType_Paramsflag) CALL MatrixTypes_Declare_ValidParams()

  !Validate against the reqParams and OptParams
  validParams=Params
  CALL validParams%validate(DenseSquareMatrixType_reqParams)

  ! Pull Data From Parameter List
  CALL validParams%get('MatrixType->n',n)
  CALL validParams%get('MatrixType->isSym',isSym)
  CALL validParams%clear()

  IF(.NOT. matrix%isInit) THEN
    IF(n < 1) THEN
      CALL eMatrixType%raiseError('Incorrect input to '// &
          modName//'::'//myName//' - Number of rows (n) must be '// &
          'greater than 1!')
    ELSE
      matrix%isInit=.TRUE.
      matrix%n=n
      IF(isSym) THEN
        matrix%isSymmetric=.TRUE.
980
      ELSE
Graham, Aaron's avatar
Graham, Aaron committed
981
        matrix%isSymmetric=.FALSE.
982
      ENDIF
Graham, Aaron's avatar
Graham, Aaron committed
983
984
985
986
987
988
989
      CALL dmallocA(matrix%a,n,n)
    ENDIF
  ELSE
    CALL eMatrixType%raiseError('Incorrect call to '// &
        modName//'::'//myName//' - MatrixType already initialized')
  ENDIF
ENDSUBROUTINE init_DenseSquareMatrixParam
990
991
992
993
994
!
!-------------------------------------------------------------------------------
!> @brief Clears the sparse matrix
!> @param matrix the matrix type to act on
!>
Graham, Aaron's avatar
Graham, Aaron committed
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
SUBROUTINE clear_SparseMatrixType(matrix)
  CLASS(SparseMatrixType),INTENT(INOUT) :: matrix
  matrix%isInit=.FALSE.
  matrix%n=0
  matrix%jCount=0
  matrix%iPrev=0
  matrix%jPrev=0
  matrix%nnz=0
  IF(ALLOCATED(matrix%ia)) CALL demallocA(matrix%ia)
  IF(ALLOCATED(matrix%ja)) CALL demallocA(matrix%ja)
  IF(ALLOCATED(matrix%a)) CALL demallocA(matrix%a)
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
ENDSUBROUTINE clear_SparseMatrixType
1008
1009
1010
1011
1012
!
!-------------------------------------------------------------------------------
!> @brief Clears the dense square matrix
!> @param matrix the matrix type to act on
!>
Graham, Aaron's avatar
Graham, Aaron committed
1013
1014
1015
1016
1017
1018
1019
1020
SUBROUTINE clear_DenseSquareMatrixType(matrix)
  CLASS(DenseSquareMatrixType),INTENT(INOUT) :: matrix
  matrix%isInit=.FALSE.
  matrix%n=0
  matrix%isSymmetric=.FALSE.
  IF(ALLOCATED(matrix%a)) CALL demallocA(matrix%a)
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
ENDSUBROUTINE clear_DenseSquareMatrixType
1021
1022
1023
1024
1025
!
!-------------------------------------------------------------------------------
!> @brief Clears the tri-diagonal matrix
!> @param matrix the matrix type to act on
!>
Graham, Aaron's avatar
Graham, Aaron committed
1026
1027
1028
1029
1030
1031
1032
1033
SUBROUTINE clear_TriDiagMatrixType(matrix)
  CLASS(TriDiagMatrixType),INTENT(INOUT) :: matrix
  matrix%isInit=.FALSE.
  matrix%n=0
  matrix%isSymmetric=.FALSE.
  IF(ALLOCATED(matrix%a)) CALL demallocA(matrix%a)
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
 ENDSUBROUTINE clear_TriDiagMatrixType
1034
1035
!
!-------------------------------------------------------------------------------
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
!> @brief Clears the banded matrix
!> @param matrix the matrix type to act on
!>
SUBROUTINE clear_BandedMatrixType(matrix)
  CLASS(BandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK) :: i

  matrix%isInit=.FALSE.
  matrix%isAssembled=.FALSE.
  matrix%n=0
  matrix%m=0
  matrix%counter=0
1048
  IF(ALLOCATED(matrix%bands)) THEN
1049
    DO i=1,SIZE(matrix%bandIdx)
1050
      IF(ALLOCATED(matrix%bands(i)%elem)) DEALLOCATE(matrix%bands(i)%elem)
1051
1052
1053
    ENDDO
    DEALLOCATE(matrix%bands)
  ENDIF
1054
1055
1056
1057
  IF(ALLOCATED(matrix%bandIdx)) DEALLOCATE(matrix%bandIdx)
  IF(ALLOCATED(matrix%iTmp)) DEALLOCATE(matrix%iTmp)
  IF(ALLOCATED(matrix%jTmp)) DEALLOCATE(matrix%jTmp)
  IF(ALLOCATED(matrix%elemTmp)) DEALLOCATE(matrix%elemTmp)
1058
  matrix%nnz=0
1059
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
 ENDSUBROUTINE clear_BandedMatrixType
!
!-------------------------------------------------------------------------------
!> @brief Clears the distributed banded matrix
!> @param matrix the matrix type to act on
!>
SUBROUTINE clear_DistributedBandedMatrixType(matrix)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: matrix
#ifdef HAVE_MPI
  INTEGER(SIK) :: i

  matrix%isInit=.FALSE.
  matrix%isCreated=.FALSE.
  matrix%isAssembled=.FALSE.
  matrix%comm=MPI_COMM_NULL
  matrix%n=0
  matrix%m=0
  matrix%nnz=0
  matrix%nLocal=0
1079
  IF(ALLOCATED(matrix%chunks)) THEN
1080
1081
1082
1083
1084
    DO i=1,SIZE(matrix%chunks)
      CALL matrix%chunks(i)%clear()
    ENDDO
    DEALLOCATE(matrix%chunks)
  ENDIF
1085
1086
1087
  IF(ALLOCATED(matrix%incIdxMap)) DEALLOCATE(matrix%incIdxMap)
  IF(ALLOCATED(matrix%incIdxStt)) DEALLOCATE(matrix%incIdxStt)
  IF(ALLOCATED(matrix%incIdxStp)) DEALLOCATE(matrix%incIdxStp)
1088
1089
1090
1091
1092
1093
  IF(ALLOCATED(matrix%iOffsets)) DEALLOCATE(matrix%iOffsets)
  IF(ALLOCATED(matrix%jOffsets)) DEALLOCATE(matrix%jOffsets)
  IF(ALLOCATED(matrix%iTmp)) DEALLOCATE(matrix%iTmp)
  IF(ALLOCATED(matrix%jTmp)) DEALLOCATE(matrix%jTmp)
  IF(ALLOCATED(matrix%elemTmp)) DEALLOCATE(matrix%elemTmp)
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
#endif
 ENDSUBROUTINE clear_DistributedBandedMatrixType

!
!-------------------------------------------------------------------------------
!> @brief Clears the distributed block banded matrix
!> @param matrix the matrix type to act on
!>
SUBROUTINE clear_DistributedBlockBandedMatrixType(matrix)
  CLASS(DistributedBlockBandedMatrixType),INTENT(INOUT) :: matrix
#ifdef HAVE_MPI
  INTEGER(SIK) :: i
  matrix%blockMask=.FALSE.
  matrix%blockSize=0
  matrix%nLocalBlocks=0
  matrix%blockOffset=0
  matrix%isInit=.FALSE.
  matrix%isCreated=.FALSE.
  matrix%isAssembled=.FALSE.

  DO i=1,SIZE(matrix%blocks)
    CALL matrix%blocks(i)%clear()
  ENDDO
  DEALLOCATE(matrix%blocks)

  CALL matrix%DistributedBandedMatrixType%clear()
#endif
 ENDSUBROUTINE clear_DistributedBlockBandedMatrixType
!
!-------------------------------------------------------------------------------
1124
1125
1126
!> @brief Clears the dense rectangular matrix
!> @param matrix the matrix type to act on
!>
Graham, Aaron's avatar
Graham, Aaron committed
1127
1128
1129
1130
1131
1132
1133
1134
SUBROUTINE clear_DenseRectMatrixType(matrix)
  CLASS(DenseRectMatrixType),INTENT(INOUT) :: matrix
  matrix%isInit=.FALSE.
  matrix%n=0
  matrix%m=0
  IF(ALLOCATED(matrix%a)) CALL demallocA(matrix%a)
  IF(MatrixType_Paramsflag) CALL MatrixTypes_Clear_ValidParams()
ENDSUBROUTINE clear_DenseRectMatrixType
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
!
!-------------------------------------------------------------------------------
!> @brief Sets the values in the sparse matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
!> This routine sets the values of the sparse matrix.  It can only be used
!> If setShape has previously been applied to the same sparse matrix.
!>
Graham, Aaron's avatar
Graham, Aaron committed
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
SUBROUTINE set_SparseMatrixType(matrix,i,j,setval)
  CLASS(SparseMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
  INTEGER(SIK) :: ja_index
  LOGICAL(SBK) :: found_ja

  REQUIRE(matrix%isInit)
  REQUIRE(matrix%jCount>0)
  REQUIRE(j>0)
  REQUIRE(i>0)
  REQUIRE(i+1<=size(matrix%ia))
  REQUIRE(i<=matrix%n)

  !currently written assuming no all-zero rows.
  !pretty safe assumption.
  found_ja=.FALSE.
  DO ja_index=matrix%ia(i),matrix%ia(i+1)-1
    IF(matrix%ja(ja_index) == j) THEN
      found_ja=.TRUE.
      EXIT
    ENDIF
  ENDDO
  IF(found_ja) matrix%a(ja_index)=setval
ENDSUBROUTINE set_SparseMatrixtype
1172
1173
!
!-------------------------------------------------------------------------------
1174
1175
1176
1177
1178
1179
1180
!> @brief Sets an entire row of values in the sparse matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j a list of j locations in the matrix
!> @param setval a list of values to set at the corresponding j locations
!>
!> This routine sets the values of an entire row in the sparse matrix.  It can only be used
Graham, Aaron's avatar
Graham, Aaron committed
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
!> If setShape has previously been applied to the same sparse matrix.
!>
SUBROUTINE setRow_SparseMatrixType(matrix,i,j,setval)
  CLASS(SparseMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j(:)
  REAL(SRK),INTENT(IN) :: setval(:)
  INTEGER(SIK) :: ja_index,ja
  LOGICAL(SBK) :: found_ja

  REQUIRE(matrix%isInit)
  REQUIRE(((matrix%jCount > 0).AND.(i <= matrix%n)) .AND. ((ALL(j > 0)) .AND. (i > 0)))
  REQUIRE(SIZE(j)==SIZE(setval))
  !currently written assuming no all-zero rows.
  !pretty safe assumption.
  DO ja = 1, SIZE(j)
    found_ja=.FALSE.
    DO ja_index=matrix%ia(i),matrix%ia(i+1)-1
      IF(matrix%ja(ja_index) == j(ja)) THEN
        found_ja=.TRUE.
        EXIT
      ENDIF
    ENDDO
    REQUIRE(found_ja)
    matrix%a(ja_index)=setval(ja)
  ENDDO
ENDSUBROUTINE setRow_SparseMatrixtype
1208
1209
1210

!
!-------------------------------------------------------------------------------
1211
1212
1213
1214
1215
1216
!> @brief Sets the values in the dense square matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
Graham, Aaron's avatar
Graham, Aaron committed
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
SUBROUTINE set_DenseSquareMatrixType(matrix,i,j,setval)
  CLASS(DenseSquareMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
  IF(matrix%isInit) THEN
    IF(((j <= matrix%n) .AND. (i <= matrix%n)) &
        .AND. ((j > 0) .AND. (i > 0))) THEN
      matrix%a(i,j)=setval
      IF(matrix%isSymmetric) matrix%a(j,i)=setval
    ENDIF
  ENDIF
ENDSUBROUTINE set_DenseSquareMatrixType
1230
1231
1232
1233
1234
1235
1236
1237
!
!-------------------------------------------------------------------------------
!> @brief Sets the values in the tridiagonal matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
Graham, Aaron's avatar
Graham, Aaron committed
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
SUBROUTINE set_TriDiagMatrixType(matrix,i,j,setval)
  CLASS(TriDiagMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
  IF(matrix%isInit) THEN
    IF(((j <= matrix%n) .AND. (i <= matrix%n)) &
        .AND. (i>=1) .AND. (j >= 1)) THEN
      !based on i,j, put in correct location
      IF((j == (i-1)).AND. (i > 1)) THEN !sub-diag
        matrix%a(1,i)=setval
        IF(matrix%isSymmetric) matrix%a(3,j)=setval
      ELSEIF((j == (i+1)) .AND. (i < matrix%n)) THEN !super-diag
        matrix%a(3,i)=setval
        IF(matrix%isSymmetric) matrix%a(1,j)=setval
      ELSEIF(i == j) THEN
        matrix%a(2,i)=setval
1255
      ENDIF
Graham, Aaron's avatar
Graham, Aaron committed
1256
1257
1258
    ENDIF
  ENDIF
ENDSUBROUTINE set_TriDiagMatrixType
1259
1260
!
!-------------------------------------------------------------------------------
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
!> @brief Sets the values in the distributed banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE set_DistributedBandedMatrixType(matrix,i,j,setval)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
#ifdef HAVE_MPI
  CHARACTER(LEN=*),PARAMETER :: myName='set_DistributedBandedMatrixType'
  INTEGER(SIK) :: rank,nproc,mpierr,k
  LOGICAL(SBK) :: flag
  CHARACTER(LEN=64) :: msg

  flag=.FALSE.
  REQUIRE(matrix%isInit)
1280
  IF(.NOT. matrix%isAssembled) THEN
1281
1282
    ! If the matrix is not yet assembled, we are adding to the temporary
    ! containers i/j/elemTmp
1283
    IF((j <= matrix%m) .AND. (i <= matrix%n) .AND. (i >= 1) .AND. (j >= 1)) THEN
1284
1285
1286
      CALL MPI_Comm_rank(matrix%comm,rank,mpierr)
      CALL MPI_Comm_size(matrix%comm,nproc,mpierr)

1287
      IF(j > matrix%jOffsets(rank+1) .AND. j <= matrix%jOffsets(rank+2)) THEN
1288
        matrix%nLocal=matrix%nLocal+1
1289
        IF(matrix%nLocal > 2*matrix%nnz/nproc) THEN
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
          CALL eMatrixType%raiseError('Matrix Element imbalance '// &
              modName//'::'//myName//' - Check nnz or use a different type')
        ENDIF
        matrix%iTmp(matrix%nLocal)=i
        matrix%jTmp(matrix%nLocal)=j
        matrix%elemTmp(matrix%nLocal)=setval
        flag=.TRUE.
      ENDIF
    ENDIF
  ELSE
    ! If the matrix is assmbled then we find the corresponding chunk and
    ! call set there
    CALL MPI_Comm_rank(matrix%comm,rank,mpierr)
1303
    IF(j > matrix%jOffsets(rank+1) .AND. j <= matrix%jOffsets(rank+2)) THEN
1304
      DO k=1,SIZE(matrix%jOffsets)-1
1305
1306
        IF(i > matrix%iOffsets(k) .AND. i <= matrix%iOffsets(k+1)) THEN
          IF(matrix%chunks(k)%isInit) THEN
1307
1308
1309
1310
1311
1312
1313
1314
1315
            CALL matrix%chunks(k)%set(i-matrix%iOffsets(k), &
                j-matrix%jOffsets(rank+1),setval)
            flag=.TRUE.
            EXIT
          ENDIF
        ENDIF
      ENDDO
    ENDIF
  ENDIF
1316
  IF(.NOT. flag) THEN
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
    WRITE(msg,'(I2,A,I9,A,I9)') rank," at ",i,',',j
    CALL eMatrixType%raiseWarning('Invalid matrix set in '// &
        modName//'::'//myName//' for rank '//msg)
  ENDIF
#endif
ENDSUBROUTINE set_DistributedBandedMatrixType

!
!-------------------------------------------------------------------------------
!> @brief Sets the values in the distributed block banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE set_DistributedBlockBandedMatrixType(matrix,i,j,setval)
  CLASS(DistributedBlockBandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
#ifdef HAVE_MPI
  CHARACTER(LEN=*),PARAMETER :: myName='set_DistributedBlockBandedMatrixType'
  INTEGER(SIK) :: rank,mpierr,offset
  LOGICAL(SBK) :: flag
  CHARACTER(LEN=64) :: msg
  REQUIRE(matrix%isInit)

  flag=.FALSE.
  CALL MPI_Comm_Rank(matrix%comm,rank,mpierr)
  ! Check if element is in a block; if not call parent
1347
  IF((i-1)/matrix%blockSize /= (j-1)/matrix%blockSize) THEN
1348
1349
1350
    CALL matrix%DistributedBandedMatrixType%set(i,j,setval)
    flag = .TRUE.
  ELSE
1351
1352
    IF((j <= matrix%m) .AND. (i <= matrix%n) .AND. (i >= 1) .AND. (j >= 1)) THEN
      IF(.NOT. matrix%blockMask) THEN
1353
        offset=((i-1)/matrix%blockSize)*matrix%blockSize
1354
        IF(j > matrix%jOffsets(rank+1) .AND. (j <= matrix%jOffsets(rank+2))) THEN
1355
1356
1357
1358
1359
1360
1361
          CALL matrix%blocks((j-1)/matrix%blockSize-matrix%blockOffset+1)%set( &
              i-offset,j-offset,setval)
          flag=.TRUE.
        ENDIF
      ENDIF
    ENDIF
  ENDIF
1362
  IF(.NOT. flag) THEN
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
    WRITE(msg,'(I2,A,I9,A,I9)') rank," at ",i,',',j
    CALL eMatrixType%raiseWarning('Invalid matrix set in '// &
        modName//'::'//myName//' for rank '//msg)
  ENDIF
#endif
ENDSUBROUTINE set_DistributedBlockBandedMatrixType

!
!-------------------------------------------------------------------------------
!> @brief Sets the values in the distributed banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE setrow_DistributedBandedMatrixType(matrix,i,j,setval)
  CHARACTER(LEN=*),PARAMETER :: myName='setrow_DistributedBandedMatrixType'
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j(:)
  REAL(SRK),INTENT(IN) :: setval(:)
  CALL eMatrixType%raiseFatalError(modName//'::'//myName// &
      ' - routine is not implemented!')
ENDSUBROUTINE setrow_DistributedBandedMatrixType
!
!-------------------------------------------------------------------------------
!> @brief Sets the values in the banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE set_BandedMatrixType(matrix,i,j,setval)
  CLASS(BandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
  INTEGER(SIK) :: bandLoc, elemIdx

  REQUIRE(matrix%isInit)
  IF((j <= matrix%m) .AND. (i <= matrix%n) .AND. (i >= 1) .AND. (j >= 1)) THEN
    IF(.NOT. matrix%isAssembled) THEN
      ! If it is not assembled, add to tmp variables
1406
      REQUIRE(matrix%counter < matrix%nnz)
1407
1408
1409
1410
1411
1412
1413
      matrix%counter=matrix%counter+1
      matrix%iTmp(matrix%counter)=i
      matrix%jTmp(matrix%counter)=j
      matrix%elemTmp(matrix%counter)=setval
    ELSE
      ! Otherwise, perform binary search to set in existing structure
      CALL binarySearch_BandedMatrixType(matrix,i,j,bandLoc,elemIdx)
1414
      IF(bandLoc >= 0 .AND. elemIdx >= 0) THEN
1415
1416
1417
1418
1419
1420
1421
        matrix%bands(bandLoc)%elem(elemIdx)=setVal
      ENDIF
    ENDIF
  ENDIF
ENDSUBROUTINE set_BandedMatrixType
!
!-------------------------------------------------------------------------------
1422
1423
1424
1425
1426
1427
!> @brief Sets the values in the dense rectangular matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
Graham, Aaron's avatar
Graham, Aaron committed
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
SUBROUTINE set_DenseRectMatrixType(matrix,i,j,setval)
  CLASS(DenseRectMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(IN) :: setval
  IF(matrix%isInit) THEN
    IF(((j <= matrix%m) .AND. (i <= matrix%n)) &
        .AND. ((j > 0) .AND. (i > 0))) matrix%a(i,j)=setval
  ENDIF
ENDSUBROUTINE set_DenseRectMatrixType
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
!
!-------------------------------------------------------------------------------
!> @brief Sets the values and the shape in the sparse matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
!> This is to be used the first time the sparse matrix values are being set.
!> This routine learns the shape of the CSR format matrix.
!> The matrix must be supplied in row-major order, any entries not in this form
!> will be ignored.
!>
Graham, Aaron's avatar
Graham, Aaron committed
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
SUBROUTINE set_shape_SparseMatrixType(matrix,i,j,setval)
  CLASS(SparseMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),OPTIONAL,INTENT(IN) :: setval

  REQUIRE(matrix%isInit)
  REQUIRE(i<=matrix%n) !row within bounds
  REQUIRE(j>0) !positive column
  REQUIRE(i>0) !positive row
  REQUIRE(matrix%jCount<matrix%nnz) !valid number of non-zero entries
  !check for increasing column order (or new row)
  REQUIRE((matrix%iPrev==i .AND. matrix%jPrev<j) .OR. (matrix%iPrev<i))

  IF(matrix%ia(i) == 0) matrix%ia(i)=matrix%jCount+1
  matrix%iPrev=i
  matrix%jPrev=j
  matrix%jCount=matrix%jCount+1
  IF(PRESENT(setval)) matrix%a(matrix%jCount)=setval
  matrix%ja(matrix%jCount)=j

ENDSUBROUTINE set_shape_SparseMatrixType
1473
1474
1475
1476
1477
1478
1479
1480
!
!-------------------------------------------------------------------------------
!> @brief Gets the values in the tridiagonal matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
Graham, Aaron's avatar
Graham, Aaron committed
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
SUBROUTINE get_TriDiagMatrixType(matrix,i,j,getval)
  CLASS(TriDiagMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(INOUT) :: getval
  IF(matrix%isInit) THEN
    IF(((j <= matrix%n) .AND. (i <= matrix%n)) &
        .AND. (i>=1) .AND. (j >= 1)) THEN
      !based on i,j, pull from correct location
      IF((j == (i-1)).AND. (i > 1)) THEN !sub-diag
        getval=matrix%a(1,i)
      ELSEIF((j == (i+1)) .AND. (i < matrix%n)) THEN !super-diag
        getval=matrix%a(3,i)
      ELSEIF(i == j) THEN
        getval=matrix%a(2,i)
1496
      ENDIF
Graham, Aaron's avatar
Graham, Aaron committed
1497
1498
1499
1500
1501
    ELSE
      getval=-1051._SRK
    ENDIF
  ENDIF
ENDSUBROUTINE get_TriDiagMatrixType
1502
1503
!
!-------------------------------------------------------------------------------
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
!> @brief Gets the values in the distributed banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE get_DistributedBandedMatrixType(matrix,i,j,getval)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(INOUT) :: getval
#ifdef HAVE_MPI
  INTEGER(SIK) :: k,mpierr,rank

  REQUIRE(matrix%isInit)
  REQUIRE(matrix%isAssembled)
  getval=0.0_SRK
  CALL MPI_Comm_rank(matrix%comm,rank,mpierr)
1522
  IF(j > matrix%jOffsets(rank+1) .AND. j <= matrix%jOffsets(rank+2)) THEN
1523
    DO k=1,SIZE(matrix%jOffsets)-1
1524
1525
      IF(i > matrix%iOffsets(k) .AND. i <= matrix%iOffsets(k+1)) THEN
        IF(matrix%chunks(k)%isInit) THEN
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
          CALL matrix%chunks(k)%get(i-matrix%iOffsets(k), &
              j-matrix%jOffsets(rank+1),getval)
          EXIT
        ENDIF
      ENDIF
    ENDDO
  ENDIF
#endif
ENDSUBROUTINE get_DistributedBandedMatrixType

!
!-------------------------------------------------------------------------------
!> @brief Gets the values in the distributed block banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE get_DistributedBlockBandedMatrixType(matrix,i,j,getval)
  CLASS(DistributedBlockBandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(INOUT) :: getval
#ifdef HAVE_MPI
  INTEGER(SIK) :: rank,mpierr,offset

  REQUIRE(matrix%isInit)
  REQUIRE(matrix%isAssembled)
  CALL MPI_Comm_Rank(matrix%comm,rank,mpierr)
1555
  IF((i-1)/matrix%blockSize /= (j-1)/matrix%blockSize) THEN
1556
    CALL matrix%DistributedBandedMatrixType%get(i,j,getVal)
1557
  ELSEIF(.NOT. matrix%blockMask) THEN
1558
    offset=((i-1)/matrix%blockSize)*matrix%blockSize
1559
    IF(j > matrix%jOffsets(rank+1) .AND. j <= matrix%jOffsets(rank+2)) THEN
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
      getval=matrix%blocks((i-1)/matrix%blockSize-matrix%blockOffset+1)%A( &
          i-offset,j-offset)
    ENDIF
  ELSE
    getval=0.0_SRK
  ENDIF
#endif
ENDSUBROUTINE get_DistributedBlockBandedMatrixType

!
!-------------------------------------------------------------------------------
!> @brief Gets the values in the banded matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
SUBROUTINE get_BandedMatrixType(matrix,i,j,getval)
  CLASS(BandedMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(INOUT) :: getval
  INTEGER(SIK) :: bandLoc,elemIdx

  REQUIRE(matrix%isInit)
  REQUIRE(matrix%isAssembled)
  CALL binarySearch_BandedMatrixType(matrix,i,j,bandLoc,elemIdx)
  getval=0.0_SRK
1588
  IF(bandLoc >= 0 .AND. elemIdx >= 0) THEN
1589
1590
1591
1592
1593
1594
    getval=matrix%bands(bandLoc)%elem(elemIdx)
  ENDIF

ENDSUBROUTINE get_BandedMatrixType
!
!-------------------------------------------------------------------------------
1595
1596
1597
1598
1599
1600
!> @brief Gets the values in the dense rectangular matrix
!> @param matrix the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!> @param setval the value to be set
!>
Graham, Aaron's avatar
Graham, Aaron committed
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
SUBROUTINE get_DenseRectMatrixType(matrix,i,j,getval)
  CLASS(DenseRectMatrixType),INTENT(INOUT) :: matrix
  INTEGER(SIK),INTENT(IN) :: i
  INTEGER(SIK),INTENT(IN) :: j
  REAL(SRK),INTENT(INOUT) :: getval
  IF(matrix%isInit) THEN
    IF(((j <= matrix%m) .AND. (i <= matrix%n)) .AND. ((j > 0) .AND. (i > 0))) THEN
      getval=matrix%a(i,j)
    ELSE
      getval=-1051._SRK
    ENDIF
  ENDIF
ENDSUBROUTINE get_DenseRectMatrixType
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
!
!-------------------------------------------------------------------------------
!> @brief Gets the values in the Dense Square matrix
!> @param declare the matrix type to act on
!> @param i the ith location in the matrix
!> @param j the jth location in the matrix
!>
!> This routine gets the values of the sparse matrix.  If the (i,j) location is
!> out of bounds, then -1051.0 (an arbitrarily chosen key) is returned.