Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module matrix_operations
5 :
6 : implicit none
7 :
8 :
9 : public :: symm_covar_matrix_2_corr_matrix, Cholesky_factor, &
10 : row_mult_lower_tri_matrix, print_lower_triangular_matrix, &
11 : get_lower_triangular_matrix, set_lower_triangular_matrix, &
12 : mirror_lower_triangular_matrix
13 :
14 : private :: Symm_matrix_eigenvalues
15 :
16 : private ! Default scope
17 :
18 : contains
19 :
20 : !-----------------------------------------------------------------------
21 0 : subroutine symm_covar_matrix_2_corr_matrix( ndim, covar, &
22 0 : corr )
23 :
24 : ! Description:
25 : ! Convert a matrix of covariances in to a matrix of correlations.
26 : ! This only does the computation the lower triangular portion of the
27 : ! matrix.
28 : ! References:
29 : ! None
30 : !-----------------------------------------------------------------------
31 :
32 : use clubb_precision, only: &
33 : core_rknd ! double precision
34 :
35 : implicit none
36 :
37 : ! External
38 : intrinsic :: sqrt
39 :
40 : ! Input Variables
41 : integer, intent(in) :: ndim
42 :
43 : real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
44 : covar ! Covariance Matrix [units vary]
45 :
46 : ! Output Variables
47 : real( kind = core_rknd ), dimension(ndim,ndim), intent(out) :: &
48 : corr ! Correlation Matrix [-]
49 :
50 : ! Local Variables
51 : integer :: i, j
52 :
53 : ! ---- Begin Code ----
54 :
55 0 : corr = 0._core_rknd ! Initialize to 0
56 :
57 0 : do i = 1, ndim
58 0 : do j = 1, i
59 0 : corr(i,j) = covar(i,j) / sqrt( covar(i,i) * covar(j,j) )
60 : end do
61 : end do
62 :
63 0 : return
64 : end subroutine symm_covar_matrix_2_corr_matrix
65 : !-----------------------------------------------------------------------
66 0 : subroutine row_mult_lower_tri_matrix( ndim, xvector, tmatrix_in, tmatrix_out )
67 :
68 : ! Description:
69 : ! Do a row-wise multiply of the elements of a lower triangular matrix.
70 : ! References:
71 : ! None
72 : !-----------------------------------------------------------------------
73 :
74 : use clubb_precision, only: &
75 : core_rknd ! double precision
76 :
77 : implicit none
78 :
79 :
80 : ! Input Variables
81 : integer, intent(in) :: ndim
82 :
83 : real( kind = core_rknd ), dimension(ndim), intent(in) :: &
84 : xvector ! Factors to be multiplied across a row [units vary]
85 :
86 : ! Input Variables
87 : real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
88 : tmatrix_in ! nxn matrix (usually a correlation matrix) [units vary]
89 :
90 : ! Output Variables
91 : real( kind = core_rknd ), dimension(ndim,ndim), intent(inout) :: &
92 : tmatrix_out ! nxn matrix (usually a covariance matrix) [units vary]
93 :
94 : ! Local Variables
95 : integer :: i, j
96 :
97 : ! ---- Begin Code ----
98 :
99 0 : do i = 1, ndim
100 0 : do j = 1, i
101 0 : tmatrix_out(i,j) = tmatrix_in(i,j) * xvector(i)
102 : end do
103 : end do
104 :
105 0 : return
106 : end subroutine row_mult_lower_tri_matrix
107 :
108 : !-------------------------------------------------------------------------------
109 0 : subroutine Cholesky_factor( ndim, a_input, a_scaling, a_Cholesky, l_scaled )
110 : ! Description:
111 : ! Create a Cholesky factorization of a_input.
112 : ! If the factorization fails we use a modified a_input matrix and attempt
113 : ! to factorize again.
114 : !
115 : ! References:
116 : ! <http://www.netlib.org/lapack/explore-html/a00868.html> dpotrf
117 : ! <http://www.netlib.org/lapack/explore-html/a00860.html> dpoequ
118 : ! <http://www.netlib.org/lapack/explore-html/a00753.html> dlaqsy
119 : !-------------------------------------------------------------------------------
120 : use error_code, only: &
121 : clubb_at_least_debug_level ! Procedure
122 :
123 : use constants_clubb, only: &
124 : fstderr ! Constant
125 :
126 : use clubb_precision, only: &
127 : core_rknd
128 :
129 : use lapack_interfaces, only: &
130 : lapack_potrf, & ! Procedures
131 : lapack_poequ, &
132 : lapack_laqsy
133 :
134 : implicit none
135 :
136 : ! Constant Parameters
137 : integer, parameter :: itermax = 10 ! Max iterations of the modified method
138 :
139 : real( kind = core_rknd), parameter :: d_coef = 0.1_core_rknd
140 : ! Coefficient applied if the decomposition doesn't work
141 :
142 : ! Input Variables
143 : integer, intent(in) :: ndim
144 :
145 : real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: a_input
146 :
147 : ! Output Variables
148 : real( kind = core_rknd ), dimension(ndim), intent(out) :: a_scaling
149 :
150 : real( kind = core_rknd ), dimension(ndim,ndim), intent(out) :: a_Cholesky
151 :
152 : logical, intent(out) :: l_scaled
153 :
154 : ! Local Variables
155 0 : real( kind = core_rknd ), dimension(ndim) :: a_eigenvalues
156 0 : real( kind = core_rknd ), dimension(ndim,ndim) :: a_corr, a_scaled
157 :
158 : real( kind = core_rknd ) :: tau, d_smallest
159 :
160 : real( kind = core_rknd ) :: amax, scond
161 : integer :: info
162 : integer :: i, j, iter
163 :
164 : character :: equed
165 :
166 : ! ---- Begin code ----
167 :
168 0 : a_scaled = a_input ! Copy input array into output array
169 :
170 : ! do i = 1, n
171 : ! do j = 1, n
172 : ! write(6,'(e10.3)',advance='no') a(i,j)
173 : ! end do
174 : ! write(6,*) ""
175 : ! end do
176 : ! pause
177 :
178 0 : equed = 'N'
179 :
180 : ! Compute scaling for a_input, using Lapack routine spoequ for single precision,
181 : ! or dpoequ for double precision
182 0 : call lapack_poequ( ndim, a_input, ndim, a_scaling, scond, amax, info )
183 :
184 0 : if ( info == 0 ) then
185 : ! Apply scaling to a_input, using Lapack routine slaqsy for single precision,
186 : ! or dlaqsy for double precision
187 0 : call lapack_laqsy( 'Lower', ndim, a_scaled, ndim, a_scaling, scond, amax, equed )
188 : end if
189 :
190 : ! Determine if scaling was necessary
191 0 : if ( equed == 'Y' ) then
192 0 : l_scaled = .true.
193 0 : a_Cholesky = a_scaled
194 : else
195 0 : l_scaled = .false.
196 0 : a_Cholesky = a_input
197 : end if
198 :
199 0 : do iter = 1, itermax
200 :
201 : ! Lapack Cholesky factorization, spotrf for single or dpotrf for double precision
202 0 : call lapack_potrf( 'Lower', ndim, a_Cholesky, ndim, info )
203 :
204 0 : select case( info )
205 : case( :-1 )
206 : write(fstderr,*) "Cholesky_factor " // &
207 0 : " illegal value for argument ", -info
208 0 : error stop
209 : case( 0 )
210 : ! Success!
211 0 : if ( clubb_at_least_debug_level( 1 ) .and. iter > 1 ) then
212 0 : write(fstderr,*) "a_factored (worked)="
213 0 : do i = 1, ndim
214 0 : do j = 1, i
215 0 : write(fstderr,'(g10.3)',advance='no') a_Cholesky(i,j)
216 : end do
217 0 : write(fstderr,*) ""
218 : end do
219 : end if
220 0 : exit
221 : case( 1: )
222 0 : if ( clubb_at_least_debug_level( 1 ) ) then
223 : ! This shouldn't happen now that the s and t Mellor(chi/eta) elements have been
224 : ! modified to never be perfectly correlated, but it's here just in case.
225 : ! -dschanen 10 Sept 2010
226 0 : write(fstderr,*) "Cholesky_factor: leading minor of order ", &
227 0 : info, " is not positive definite."
228 0 : write(fstderr,*) "factorization failed."
229 0 : write(fstderr,*) "a_input="
230 0 : do i = 1, ndim
231 0 : do j = 1, i
232 0 : write(fstderr,'(g10.3)',advance='no') a_input(i,j)
233 : end do
234 0 : write(fstderr,*) ""
235 : end do
236 0 : write(fstderr,*) "a_Cholesky="
237 0 : do i = 1, ndim
238 0 : do j = 1, i
239 0 : write(fstderr,'(g10.3)',advance='no') a_Cholesky(i,j)
240 : end do
241 0 : write(fstderr,*) ""
242 : end do
243 : end if
244 :
245 0 : if ( clubb_at_least_debug_level( 2 ) ) then
246 : call Symm_matrix_eigenvalues( ndim, a_input, & ! intent(in)
247 0 : a_eigenvalues ) ! intent(out)
248 0 : write(fstderr,*) "a_eigenvalues="
249 0 : do i = 1, ndim
250 0 : write(fstderr,'(g10.3)',advance='no') a_eigenvalues(i)
251 : end do
252 0 : write(fstderr,*) ""
253 :
254 : call symm_covar_matrix_2_corr_matrix( ndim, a_input, & ! intent(in)
255 0 : a_corr ) ! intent(out)
256 0 : write(fstderr,*) "a_correlations="
257 0 : do i = 1, ndim
258 0 : do j = 1, i
259 0 : write(fstderr,'(g10.3)',advance='no') a_corr(i,j)
260 : end do
261 0 : write(fstderr,*) ""
262 : end do
263 : end if
264 :
265 0 : if ( iter == itermax ) then
266 0 : write(fstderr,*) "iteration =", iter, "itermax =", itermax
267 0 : write(fstderr,*) "Fatal error in Cholesky_factor"
268 0 : else if ( clubb_at_least_debug_level( 1 ) ) then
269 : ! Adding a STOP statement to prevent this problem from slipping under
270 : ! the rug.
271 0 : write(fstderr,*) "Fatal error in Cholesky_factor"
272 0 : write(fstderr,*) "Attempting to modify matrix to allow factorization."
273 : end if
274 :
275 0 : if ( l_scaled ) then
276 0 : a_Cholesky = a_scaled
277 : else
278 0 : a_Cholesky = a_input
279 : end if
280 : ! The number used for tau here is case specific to the Sigma covariance
281 : ! matrix in the latin hypercube code and is not at all general.
282 : ! Tau should be a number that is small relative to the other diagonal
283 : ! elements of the matrix to have keep the error caused by modifying 'a' low.
284 : ! -dschanen 30 Aug 2010
285 0 : d_smallest = a_Cholesky(1,1)
286 0 : do i = 2, ndim
287 0 : if ( d_smallest > a_Cholesky(i,i) ) d_smallest = a_Cholesky(i,i)
288 : end do
289 : ! Use the smallest element * d_coef * iteration
290 0 : tau = d_smallest * d_coef * real( iter, kind=core_rknd )
291 :
292 : ! print *, "tau =", tau, "d_smallest = ", d_smallest
293 :
294 0 : do i = 1, ndim
295 0 : do j = 1, ndim
296 0 : if ( i == j ) then
297 0 : a_Cholesky(i,j) = a_Cholesky(i,j) + tau ! Add tau to the diagonal
298 : else
299 0 : a_Cholesky(i,j) = a_Cholesky(i,j)
300 : end if
301 : end do
302 : end do
303 :
304 0 : if ( clubb_at_least_debug_level( 2 ) ) then
305 : call Symm_matrix_eigenvalues( ndim, a_Cholesky, & ! intent(in)
306 0 : a_eigenvalues ) ! intent(out)
307 0 : write(fstderr,*) "a_modified eigenvalues="
308 0 : do i = 1, ndim
309 0 : write(fstderr,'(e10.3)',advance='no') a_eigenvalues(i)
310 : end do
311 0 : write(fstderr,*) ""
312 : end if
313 :
314 : end select ! info
315 : end do ! 1..itermax
316 :
317 0 : return
318 : end subroutine Cholesky_factor
319 :
320 : !----------------------------------------------------------------------
321 0 : subroutine Symm_matrix_eigenvalues( ndim, a_input, &
322 0 : a_eigenvalues )
323 :
324 : ! Description:
325 : ! Computes the eigevalues of a_input
326 : !
327 : ! References:
328 : ! None
329 : !-----------------------------------------------------------------------
330 :
331 : use constants_clubb, only: &
332 : fstderr ! Constant
333 :
334 : use clubb_precision, only: &
335 : core_rknd ! double precision
336 :
337 : use lapack_interfaces, only: &
338 : lapack_syev ! Procedure
339 :
340 : implicit none
341 :
342 : ! Parameters
343 : integer, parameter :: &
344 : lwork = 180 ! This is the optimal value I obtained for an n of 5 -dschanen 31 Aug 2010
345 :
346 : ! Input Variables
347 : integer, intent(in) :: ndim
348 :
349 : real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: a_input
350 :
351 : ! Output Variables
352 : real( kind = core_rknd ), dimension(ndim), intent(out) :: a_eigenvalues
353 :
354 : ! Local Variables
355 0 : real( kind = core_rknd ), dimension(ndim,ndim) :: a_scratch
356 :
357 : real( kind = core_rknd ), dimension(lwork) :: work
358 :
359 : integer :: info
360 : ! integer :: i, j
361 : ! ---- Begin code ----
362 :
363 0 : a_scratch = a_input
364 :
365 : ! do i = 1, ndim
366 : ! do j = 1, ndim
367 : ! write(6,'(e10.3)',advance='no') a(i,j)
368 : ! end do
369 : ! write(6,*) ""
370 : ! end do
371 : ! pause
372 :
373 : ! Lapack routine for computing eigenvalues and, optionally, eigenvectors. ssyev for
374 : ! single precision or dsyev for double precision
375 : call lapack_syev( 'No eigenvectors', 'Lower', ndim, a_scratch, ndim, &
376 0 : a_eigenvalues, work, lwork, info )
377 :
378 0 : select case( info )
379 : case( :-1 )
380 : write(fstderr,*) "Symm_matrix_eigenvalues:" // &
381 0 : " illegal value for argument ", -info
382 : case( 0 )
383 : ! Success!
384 :
385 : case( 1: )
386 0 : write(fstderr,*) "Symm_matrix_eigenvalues: Algorithm failed to converge."
387 : end select
388 :
389 0 : return
390 : end subroutine Symm_matrix_eigenvalues
391 : !-------------------------------------------------------------------------------
392 0 : subroutine set_lower_triangular_matrix( pdf_dim, index1, index2, xpyp, &
393 0 : matrix )
394 : ! Description:
395 : ! Set a value for the lower triangular portion of a matrix.
396 : ! References:
397 : ! None
398 : !-------------------------------------------------------------------------------
399 :
400 : use clubb_precision, only: &
401 : core_rknd ! user defined precision
402 :
403 : implicit none
404 :
405 : ! External
406 : intrinsic :: max, min
407 :
408 : ! Input Variables
409 : integer, intent(in) :: &
410 : pdf_dim, & ! Number of variates
411 : index1, index2 ! Indices for 2 variates (the order doesn't matter)
412 :
413 : real( kind = core_rknd ), intent(in) :: &
414 : xpyp ! Value for the matrix (usually a correlation or covariance) [units vary]
415 :
416 : ! Input/Output Variables
417 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(inout) :: &
418 : matrix ! The lower triangular matrix
419 :
420 : integer :: i,j
421 :
422 : ! ---- Begin Code ----
423 :
424 : ! Reverse these to set the values of upper triangular matrix
425 0 : i = max( index1, index2 )
426 0 : j = min( index1, index2 )
427 :
428 0 : if( i > 0 .and. j > 0 ) then
429 0 : matrix(i,j) = xpyp
430 : end if
431 :
432 0 : return
433 : end subroutine set_lower_triangular_matrix
434 : !-------------------------------------------------------------------------------
435 :
436 : !-------------------------------------------------------------------------------
437 0 : subroutine get_lower_triangular_matrix( pdf_dim, index1, index2, matrix, &
438 : xpyp )
439 : ! Description:
440 : ! Returns a value from the lower triangular portion of a matrix.
441 : ! References:
442 : ! None
443 : !-------------------------------------------------------------------------------
444 :
445 : use clubb_precision, only: &
446 : core_rknd
447 :
448 : implicit none
449 :
450 : ! External
451 : intrinsic :: max, min
452 :
453 : ! Input Variables
454 : integer, intent(in) :: &
455 : pdf_dim, & ! Number of variates
456 : index1, index2 ! Indices for 2 variates (the order doesn't matter)
457 :
458 : ! Input/Output Variables
459 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(in) :: &
460 : matrix ! The covariance matrix
461 :
462 : real( kind = core_rknd ), intent(out) :: &
463 : xpyp ! Value from the matrix (usually a correlation or covariance) [units vary]
464 :
465 : integer :: i,j
466 :
467 : ! ---- Begin Code ----
468 :
469 : ! Reverse these to set the values of upper triangular matrix
470 0 : i = max( index1, index2 )
471 0 : j = min( index1, index2 )
472 :
473 0 : xpyp = matrix(i,j)
474 :
475 0 : return
476 : end subroutine get_lower_triangular_matrix
477 :
478 : !-----------------------------------------------------------------------
479 0 : subroutine print_lower_triangular_matrix( iunit, ndim, matrix )
480 :
481 : ! Description:
482 : ! Print the values of lower triangular matrix to a file or console.
483 :
484 : ! References:
485 : ! None
486 : !-----------------------------------------------------------------------
487 :
488 : use clubb_precision, only: &
489 : core_rknd ! Variable(s)
490 :
491 : implicit none
492 :
493 : ! Input Variables
494 : integer, intent(in) :: &
495 : iunit, & ! File I/O logical unit (usually 6 for stdout and 0 for stderr)
496 : ndim ! Dimension of the matrix
497 :
498 : real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
499 : matrix ! Lower triangular matrix [units vary]
500 :
501 : ! Local Variables
502 : integer :: i, j
503 :
504 : ! ---- Begin Code ----
505 :
506 0 : do i = 1, ndim
507 0 : do j = 1, i
508 0 : write(iunit,fmt='(g15.6)',advance='no') matrix(i,j)
509 : end do
510 0 : write(iunit,fmt=*) "" ! newline
511 : end do
512 :
513 0 : return
514 : end subroutine print_lower_triangular_matrix
515 :
516 : !-----------------------------------------------------------------------
517 0 : subroutine mirror_lower_triangular_matrix( nvars, &
518 0 : matrix )
519 :
520 : ! Description:
521 : ! Mirrors the elements of a lower triangular matrix to the upper
522 : ! triangle so that it is symmetric.
523 :
524 : ! References:
525 : ! None
526 : !-----------------------------------------------------------------------
527 :
528 : use clubb_precision, only: &
529 : core_rknd ! Constant
530 :
531 : implicit none
532 :
533 : ! Input Variables
534 : integer, intent(in) :: &
535 : nvars ! Number of variables in each dimension of square matrix
536 :
537 : ! Input/Output Variables
538 : real( kind = core_rknd ), dimension(nvars,nvars), intent(inout) :: &
539 : matrix ! Lower triangluar square matrix
540 :
541 : ! Local Variables
542 : integer :: row, col
543 :
544 : !-----------------------------------------------------------------------
545 :
546 : !----- Begin Code -----
547 :
548 0 : if ( nvars > 1 ) then
549 :
550 0 : do col=2, nvars
551 0 : do row=1, col-1
552 :
553 0 : matrix(row,col) = matrix(col,row)
554 :
555 : end do
556 : end do
557 :
558 : end if ! nvars > 1
559 :
560 0 : return
561 :
562 : end subroutine mirror_lower_triangular_matrix
563 : !-----------------------------------------------------------------------
564 :
565 : end module matrix_operations
|
