ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_matrixChol::setMatChol Interface Reference

Compute and return the lower/upper-triangle of the Cholesky factorization \(L\) of the input Symmetric/Hermitian positive-definite triangular matrix.
More...

Detailed Description

Compute and return the lower/upper-triangle of the Cholesky factorization \(L\) of the input Symmetric/Hermitian positive-definite triangular matrix.

Parameters
[in,out]mat: The input or input/output array of rank 2 of square shape (1 : ndim, 1 : ndim) of,
  1. type complex of kind any supported by the processor (e.g., CK, CK32, CK64, or CK128), or
  2. type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128).
The intent of the argument mat depends on the presence or settings of the rest of the arguments of the generic interface.
  1. On input, the specified upper/lower-diagonal subset of mat must contain the square symmetric positive-definite matrix whose Cholesky factorization is to be computed.
  2. On output,
    1. If the optional output argument chol is missing, then mat has intent(inout) and the same specified triangular subset of mat will contain the Cholesky factorization of the matrix, while the other triangle of mat remains intact.
      In this case, the blocking algorithms of BLAS/LAPACK are used for the factorization.
    2. If the optional output argument chol is present, then the input argument mat has intent(in) and remains intact on return.
      The specified triangular subset of mat will be used to construct the same (or the opposite) triangular subset of the Cholesky factorization of mat.
      The exact storage format of the output chol is determined by the input arguments subset and operation.
      In all cases, the other side of chol remains intact upon return.
      For example,
      1. setting subset = uppDia, operation = nothing will lead to reading the matrix from the upper-diagonal triangular subset of mat and writing the upper-triangular Cholesky factorization to the upper-diagonal triangular subset of chol while its upper triangle remains intact.
      2. setting subset = uppDia, operation = transHerm will lead to reading the matrix from the upper-diagonal triangular subset of mat and writing the lower-triangular Cholesky factorization to the lower-diagonal triangular subset of chol while its upper triangle remains intact.
      3. setting subset = lowDia, operation = nothing will lead to reading the matrix from the lower-diagonal triangular subset of mat and writing the corresponding Cholesky factorization to the upper-diagonal triangular subset of chol while its upper triangle remains intact.
      4. setting subset = lowDia, operation = transHerm will lead to reading the matrix from the lower-diagonal triangular subset of mat and writing the lower-triangular Cholesky factorization to the lower-diagonal triangular subset of chol while its upper triangle remains intact.
      The input mat and chol can be two opposite subsets of the same matrix.
      This format is particularly useful for efficient packing of mat with its Cholesky factorization with losing the diagonals of either matrix.
      The default algorithms with this format are the non-blocking basic (Cholesky–Crout / Cholesky–Banachiewicz) algorithms.
      These algorithms are very fast for low rank input mat but slower than the blocked LAPACK algorithms at higher matrix ranks.
[in]subset: The input scalar that can be either,
  1. the constant uppDia implying that the upper-diagonal triangular block of mat should be used for building the Cholesky factor while the lower triangular part of mat remains intact on output.
  2. the constant lowDia implying that the lower-diagonal triangular block of mat should be used for building the Cholesky factor while the upper triangular part of mat remains intact on output.
This argument is merely a convenience to differentiate the different procedure functionalities within this generic interface.
[out]info: The output scalar integer of default kind IK that is 0 if and only if the Cholesky factorization succeeds.
Otherwise, it is set to the order of the leading minor of the specified input subset of mat that is not positive definite, indicating the occurrence of an error and that the factorization could not be completed.
A non-zero value implies failure of the Cholesky factorization.
[in,out]chol: The input/output matrix of the same type, kind, and shape as the input mat containing the Cholesky factorization in its triangular subset dictated by the input arguments subset and operation.
See the documentation of the input argument mat above for more information.
(optional. It must be present if and only if the input argument operation is present and control is missing.)
[in]operation: The input scalar that can be either,
  1. the constant nothing implying that the same subset of the output chol as the input subset must be populated with the Cholesky factorization.
  2. the constant transHerm implying that the complex transpose of the computed Cholesky factorization must be written to the output chol.
(optional, default = nothing. It must be present if and only if the output argument chol is also present.)
[in]control: The input scalar that can be,
  1. the constant recursion implying that the recursive version of the algorithm must be used.
(optional. It can be present only if the input arguments chol and operation are both missing.)
[in]ndim: The input scalar integer of default kind IK representing the order of the square subset (roff : roff + ndim, coff : coff + ndim) of the input matrix of arbitrary shape (:,:).
(optional, default = size(mat, 1).)
[in]roff: The input non-negative scalar of type integer of default kind IK containing the offset from the first row of the input matrix mat, such that mat(1 + roff, 1 + roff) marks the top-left corner of the block of mat used in the algorithm.
(optional, default = 0.)
[in]coff: The input non-negative scalar of type integer of default kind IK containing the offset from the first column of the input matrix mat, such that mat(1 + roff, 1 + roff) marks the top-left corner of the block of mat used in the algorithm.
(optional, default = 0.)
[in]roffc: The input non-negative scalar of type integer of default kind IK containing the offset from the first row of the output matrix chol, such that chol(1 + roffc, 1 + roffc) marks the top-left corner of the block of chol used in the algorithm.
(optional, default = 0.)
[in]coffc: The input non-negative scalar of type integer of default kind IK containing the offset from the first column of the output matrix chol, such that chol(1 + roffc, 1 + roffc) marks the top-left corner of the block of chol used in the algorithm.
(optional, default = 0.)
[in]bdim: The input scalar integer of default kind IK representing the optimal block size to be used in the blocked algorithm.
For any input ndim <= bdim, the procedure uses the unblocked recursive algorithm, otherwise, the blocked algorithm.
(optional. default = 64. It can be present only if the input argument control is missing.)


Possible calling interfaces

! simplified interface.
call setMatChol(mat, subset, info)
call setMatChol(mat, subset, info, control) ! control = recursion
call setMatChol(mat, subset, info, control, bdim = bdim) ! control = iteration
call setMatChol(mat, subset, info, chol, operation)
! contiguous interface.
call setMatChol(mat, subset, info, ndim, roff, coff)
call setMatChol(mat, subset, info, control, ndim, roff, coff) ! control = recursion
call setMatChol(mat, subset, info, control, ndim, roff, coff, bdim = bdim) ! control = iteration
call setMatChol(mat, subset, info, chol, operation, ndim, roff, coff, roffc, coffc)
Compute and return the lower/upper-triangle of the Cholesky factorization of the input Symmetric/Her...
This module contains procedures and generic interfaces for computing the Cholesky factorization of po...
Warning
The condition 0 <= ndim must hold for the corresponding input arguments.
The condition all(0 <= [roff, coff]) must hold for the corresponding input arguments.
The condition all([roff, coff] + ndim <= shape(mat)) must hold for the corresponding input arguments.
The condition all(shape(chol) == shape(mat)) must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.
The pure procedure(s) documented herein become impure when the ParaMonte library is compiled with preprocessor macro CHECK_ENABLED=1.
By default, these procedures are pure in release build and impure in debug and testing builds.
See also
pm_matrixDet
pm_matrixInv
pm_matrixLUP
transHerm
iteration
recursion
lowDia
uppDia
BLAS/LAPACK equivalent:
The procedures under discussion combine, modernize, and extend the interface and functionalities of Version 3.11 of BLAS/LAPACK routine(s): SPOTRF2, DPOTRF2, CPOTRF2, and ZPOTRF2.


Example usage

1program example
2
3 use pm_kind, only: SK, IK, LK
4 use pm_kind, only: TKG => RKS ! all processor type kinds are supported.
5 use pm_kind, only: CKG => CKS ! all processor type kinds are supported.
6 use pm_matrixChol, only: setMatChol, lowDia, uppDia
7 use pm_matrixChol, only: iteration, recursion
8 use pm_matrixChol, only: nothing, transHerm
10 use pm_io, only: display_type
11 use pm_io, only: getFormat
12
13 implicit none
14
15 type(display_type) :: disp
16
17 character(:, SK), allocatable :: cform, gform
18 real(TKG) , parameter :: DUM = -huge(0._TKG)
19 complex(CKG) , parameter :: CMPLX_DUMM = cmplx(-huge(0._CKG), -huge(0._CKG), CKG)
20 integer(IK) :: info, ndim, roff, coff
21 cform = getFormat([CMPLX_DUMM], ed = SK_'f', signed = .true.)
22 gform = getFormat([DUM], ed = SK_'f', signed = .true.)
23
24 disp = display_type(file = "main.out.F90")
25
26 block
27 use pm_kind, only: TKG => RKS
28 real(TKG), allocatable :: mat(:,:), cholow(:,:), choupp(:,:)
29 mat = reshape( [ 1._TKG, 0._TKG, 2._TKG &
30 , 0._TKG, 4._TKG, 0._TKG &
31 , 2._TKG, 0._TKG, 8._TKG &
32 ], shape = [3,3], order = [2, 1])
33 call disp%skip
34 call disp%show("mat")
35 call disp%show( mat )
36 call setResized(choupp, shape(mat, IK))
37 call disp%show("choupp = 0")
38 choupp = 0
39 call disp%show("call setMatChol(mat, uppDia, info, choupp, nothing)")
40 call setMatChol(mat, uppDia, info, choupp, nothing)
41 call disp%show("if (info /= 0) error stop")
42 if (info /= 0) error stop
43 call disp%show("choupp")
44 call disp%show( choupp )
45 call setResized(cholow, shape(mat, IK))
46 call disp%show("cholow = 0")
47 cholow = 0
48 call disp%show("call setMatChol(mat, uppDia, info, cholow, transHerm)")
49 call setMatChol(mat, uppDia, info, cholow, transHerm)
50 call disp%show("if (info /= 0) error stop")
51 if (info /= 0) error stop
52 call disp%show("cholow")
53 call disp%show( cholow )
54 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
55 call disp%show( matmul(cholow, choupp) - mat )
56 call disp%skip
57 end block
58
59 block
60 use pm_kind, only: TKG => RKS
61 complex(TKG), allocatable :: mat(:,:), cholow(:,:), choupp(:,:)
62 mat = reshape( [ (9.0, 0.0), (3.0, 3.0), (3.0, -3.0) &
63 , (3.0, -3.0),(18.0, 0.0), (8.0, -6.0) &
64 , (3.0, 3.0), (8.0, 6.0),(43.0, 0.0) &
65 ], shape = [3,3], order = [2, 1])
66 call disp%skip
67 call disp%show("mat")
68 call disp%show( mat )
69 call setResized(choupp, shape(mat, IK))
70 call disp%show("choupp = 0")
71 choupp = 0
72 call disp%show("call setMatChol(mat, uppDia, info, choupp, nothing)")
73 call setMatChol(mat, uppDia, info, choupp, nothing)
74 call disp%show("if (info /= 0) error stop")
75 if (info /= 0) error stop
76 call disp%show("choupp")
77 call disp%show( choupp )
78 call setResized(cholow, shape(mat, IK))
79 call disp%show("cholow = 0")
80 cholow = 0
81 call disp%show("call setMatChol(mat, uppDia, info, cholow, transHerm)")
82 call setMatChol(mat, uppDia, info, cholow, transHerm)
83 call disp%show("if (info /= 0) error stop")
84 if (info /= 0) error stop
85 call disp%show("cholow")
86 call disp%show( cholow )
87 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
88 call disp%show( matmul(cholow, choupp) - mat )
89 call disp%skip
90 end block
91
92 block
93 use pm_kind, only: TKG => RKS
94 complex(TKG), allocatable :: mat(:,:), cholow(:,:), choupp(:,:)
95 mat = reshape( [ (25.0, 0.0), (-5.0, -5.0), (10.0, 5.0) &
96 , (-5.0, 5.0), (51.0, 0.0), (4.0, -6.0) &
97 , (10.0, -5.0), (4.0, 6.0), (71.0, 0.0) &
98 ], shape = [3,3], order = [2, 1])
99 call disp%skip
100 call disp%show("mat")
101 call disp%show( mat )
102 call setResized(choupp, shape(mat, IK))
103 call disp%show("choupp = 0")
104 choupp = 0
105 call disp%show("call setMatChol(mat, uppDia, info, choupp, nothing)")
106 call setMatChol(mat, uppDia, info, choupp, nothing)
107 call disp%show("if (info /= 0) error stop")
108 if (info /= 0) error stop
109 call disp%show("choupp")
110 call disp%show( choupp )
111 call setResized(cholow, shape(mat, IK))
112 call disp%show("cholow = 0")
113 cholow = 0
114 call disp%show("call setMatChol(mat, uppDia, info, cholow, transHerm)")
115 call setMatChol(mat, uppDia, info, cholow, transHerm)
116 call disp%show("if (info /= 0) error stop")
117 if (info /= 0) error stop
118 call disp%show("cholow")
119 call disp%show( cholow )
120 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
121 call disp%show( matmul(cholow, choupp) - mat )
122 call disp%skip
123 end block
124
125 block
126 use pm_kind, only: TKG => RKS
127 real(TKG), allocatable :: mat(:,:), choupp(:,:), cholow(:,:)
128 mat = reshape( [ 1._TKG, 0._TKG, 2._TKG &
129 , 0._TKG, 4._TKG, 0._TKG &
130 , 2._TKG, 0._TKG, 8._TKG &
131 ], shape = [3,3], order = [2, 1])
132 call disp%skip
133 call disp%show("mat")
134 call disp%show( mat )
135 call setResized(cholow, shape(mat, IK))
136 call disp%show("cholow = 0")
137 cholow = 0
138 call disp%show("call setMatChol(mat, lowDia, info, cholow, nothing)")
139 call setMatChol(mat, lowDia, info, cholow, nothing)
140 call disp%show("if (info /= 0) error stop")
141 if (info /= 0) error stop
142 call disp%show("cholow")
143 call disp%show( cholow )
144 call setResized(choupp, shape(mat, IK))
145 call disp%show("choupp = 0")
146 call disp%show("choupp = 0")
147 choupp = 0
148 call disp%show("call setMatChol(mat, lowDia, info, choupp, transHerm)")
149 call setMatChol(mat, lowDia, info, choupp, transHerm)
150 call disp%show("if (info /= 0) error stop")
151 if (info /= 0) error stop
152 call disp%show("choupp")
153 call disp%show( choupp )
154 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
155 call disp%show( matmul(cholow, choupp) - mat )
156 call disp%skip
157 end block
158
159 block
160 use pm_kind, only: TKG => RKS
161 complex(TKG), allocatable :: mat(:,:), choupp(:,:), cholow(:,:)
162 mat = reshape( [ (9.0, 0.0), (3.0, 3.0), (3.0, -3.0) &
163 , (3.0, -3.0),(18.0, 0.0), (8.0, -6.0) &
164 , (3.0, 3.0), (8.0, 6.0),(43.0, 0.0) &
165 ], shape = [3,3], order = [2, 1])
166 call disp%skip
167 call disp%show("mat")
168 call disp%show( mat )
169 call setResized(cholow, shape(mat, IK))
170 call disp%show("cholow = 0")
171 cholow = 0
172 call disp%show("call setMatChol(mat, lowDia, info, cholow, nothing)")
173 call setMatChol(mat, lowDia, info, cholow, nothing)
174 call disp%show("if (info /= 0) error stop")
175 if (info /= 0) error stop
176 call disp%show("cholow")
177 call disp%show( cholow )
178 call setResized(choupp, shape(mat, IK))
179 call disp%show("choupp = 0")
180 call disp%show("choupp = 0")
181 choupp = 0
182 call disp%show("call setMatChol(mat, lowDia, info, choupp, transHerm)")
183 call setMatChol(mat, lowDia, info, choupp, transHerm)
184 call disp%show("if (info /= 0) error stop")
185 if (info /= 0) error stop
186 call disp%show("choupp")
187 call disp%show( choupp )
188 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
189 call disp%show( matmul(cholow, choupp) - mat )
190 call disp%skip
191 end block
192
193 block
194 use pm_kind, only: TKG => RKS
195 complex(TKG), allocatable :: mat(:,:), choupp(:,:), cholow(:,:)
196 mat = reshape( [ (25.0, 0.0), (-5.0, -5.0), (10.0, 5.0) &
197 , (-5.0, 5.0), (51.0, 0.0), (4.0, -6.0) &
198 , (10.0, -5.0), (4.0, 6.0), (71.0, 0.0) &
199 ], shape = [3,3], order = [2, 1])
200 call disp%skip
201 call disp%show("mat")
202 call disp%show( mat )
203 call setResized(cholow, shape(mat, IK))
204 call disp%show("cholow = 0")
205 cholow = 0
206 call disp%show("call setMatChol(mat, lowDia, info, cholow, nothing)")
207 call setMatChol(mat, lowDia, info, cholow, nothing)
208 call disp%show("if (info /= 0) error stop")
209 if (info /= 0) error stop
210 call disp%show("cholow")
211 call disp%show( cholow )
212 call setResized(choupp, shape(mat, IK))
213 call disp%show("choupp = 0")
214 call disp%show("choupp = 0")
215 choupp = 0
216 call disp%show("call setMatChol(mat, lowDia, info, choupp, transHerm)")
217 call setMatChol(mat, lowDia, info, choupp, transHerm)
218 call disp%show("if (info /= 0) error stop")
219 if (info /= 0) error stop
220 call disp%show("choupp")
221 call disp%show( choupp )
222 call disp%show("matmul(cholow, choupp) - mat ! must be all zero-valued elements.")
223 call disp%show( matmul(cholow, choupp) - mat )
224 call disp%skip
225 end block
226
227 call disp%skip
228 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
229 call disp%show("!In-place and out-of-place factorization within the same matrix")
230 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
231 call disp%skip
232
233 block
234 integer(IK), parameter :: ndim = 9
235 real(TKG) :: cholref(ndim, ndim), mat(ndim, ndim), matchol(ndim, ndim + 1)
236 mat = reshape([real(TKG) ::1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
237 , 1.0, 2.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
238 , 1.0, 2.0, 3.0, DUM, DUM, DUM, DUM, DUM, DUM &
239 , 1.0, 2.0, 3.0, 4.0, DUM, DUM, DUM, DUM, DUM &
240 , 1.0, 2.0, 3.0, 4.0, 5.0, DUM, DUM, DUM, DUM &
241 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, DUM, DUM, DUM &
242 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, DUM, DUM &
243 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, DUM &
244 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 &
245 ], shape = [ndim, ndim], order = [2, 1])
246 cholref = reshape([real(TKG)::1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
247 , 1.0, 1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
248 , 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM, DUM, DUM &
249 , 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM, DUM &
250 , 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM &
251 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM &
252 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM &
253 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM &
254 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
255 ], shape = [ndim, ndim], order = [2, 1])
256 call disp%skip()
257 call disp%show("mat")
258 call disp%show( mat , format = gform )
259 matchol = DUM
260 call disp%show("matchol(:,1:ndim) = mat")
261 matchol(:,1:ndim) = mat
262 call disp%show("matchol(:,1:ndim)")
263 call disp%show( matchol(:,1:ndim) , format = gform )
264 call disp%show("call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)")
265 call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)
266 call disp%show("matchol(:,1:ndim)")
267 call disp%show( matchol(:,1:ndim) , format = gform )
268 call disp%show("matchol(:,1:ndim) - cholref")
269 call disp%show( matchol(:,1:ndim) - cholref , format = gform )
270 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
271 if (info /= 0) error stop 'Cholesky factorization failed.'
272 call disp%skip()
273 matchol = DUM
274 call disp%show("matchol(:,1:ndim) = mat")
275 matchol(:,1:ndim) = mat
276 call disp%show("matchol(:,1:ndim)")
277 call disp%show( matchol(:,1:ndim) , format = gform )
278 call disp%show("call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)")
279 call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)
280 call disp%show("matchol(:,2:ndim+1)")
281 call disp%show( matchol(:,2:ndim+1) , format = gform )
282 call disp%show("matchol(:,2:ndim+1) - transpose(cholref)")
283 call disp%show( matchol(:,2:ndim+1) - transpose(cholref) , format = gform )
284 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
285 if (info /= 0) error stop 'Cholesky factorization failed.'
286 call disp%skip()
287
288 end block
289
290 block
291 integer(IK), parameter :: ndim = 3
292 complex(TKG) :: cholref(ndim, ndim), mat(ndim, ndim), matchol(ndim, ndim + 1)
293 mat = reshape( [complex(CKG) ::(25.0, 0.0), CMPLX_DUMM, CMPLX_DUMM &
294 , (-5.0, 5.0), (51.0, 0.0), CMPLX_DUMM &
295 , (10.0, -5.0), ( 4.0, 6.0), (71.0, 0.0) &
296 ], shape = [ndim, ndim], order = [2, 1])
297 cholref = reshape( [complex(CKG)::( 5.0, 0.0), CMPLX_DUMM, CMPLX_DUMM &
298 , (-1.0, 1.0), (7.0, 0.0), CMPLX_DUMM &
299 , ( 2.0, -1.0), (1.0, 1.0), (8.0, 0.0) &
300 ], shape = [ndim, ndim], order = [2, 1])
301 call disp%skip()
302 call disp%show("mat")
303 call disp%show( mat , format = cform )
304 matchol = CMPLX_DUMM
305 call disp%show("matchol(:,1:ndim) = mat")
306 matchol(:,1:ndim) = mat
307 call disp%show("matchol(:,1:ndim)")
308 call disp%show( matchol(:,1:ndim) , format = cform )
309 call disp%show("call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)")
310 call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)
311 call disp%show("matchol(:,1:ndim)")
312 call disp%show( matchol(:,1:ndim) , format = cform )
313 call disp%show("cholref")
314 call disp%show( cholref , format = cform )
315 call disp%show("matchol(:,1:ndim) - cholref")
316 call disp%show( matchol(:,1:ndim) - cholref , format = cform )
317 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
318 if (info /= 0) error stop 'Cholesky factorization failed.'
319 call disp%skip()
320 matchol = CMPLX_DUMM
321 call disp%show("matchol(:,1:ndim) = mat")
322 matchol(:,1:ndim) = mat
323 call disp%show("matchol(:,1:ndim)")
324 call disp%show( matchol(:,1:ndim) , format = cform )
325 call disp%show("call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)")
326 call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)
327 call disp%show("matchol(:,2:ndim+1)")
328 call disp%show( matchol(:,2:ndim+1) , format = cform )
329 call disp%show("transpose(conjg(cholref))")
330 call disp%show( transpose(conjg(cholref)) , format = cform )
331 call disp%show("matchol(:,2:ndim+1) - transpose(conjg(cholref))")
332 call disp%show( matchol(:,2:ndim+1) - transpose(conjg(cholref)) , format = cform )
333 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
334 if (info /= 0) error stop 'Cholesky factorization failed.'
335 call disp%skip()
336
337 end block
338
339 call runExamples()
340 call runExamples(bdim = 2_IK)
341
342contains
343
344 subroutine runExamples(bdim)
345 integer(IK), optional :: bdim
346
347 call disp%skip
348 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
349 call disp%show("!Cholesky factorization of real positive-definite square matrix update.")
350 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
351 call disp%skip
352
353 block
354
355 real(TKG), allocatable :: cholref(:,:), cholmat(:,:)
356
357 cholmat = reshape([real(TKG) :: 1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
358 , 1.0, 2.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
359 , 1.0, 2.0, 3.0, DUM, DUM, DUM, DUM, DUM, DUM &
360 , 1.0, 2.0, 3.0, 4.0, DUM, DUM, DUM, DUM, DUM &
361 , 1.0, 2.0, 3.0, 4.0, 5.0, DUM, DUM, DUM, DUM &
362 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, DUM, DUM, DUM &
363 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, DUM, DUM &
364 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, DUM &
365 , 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 &
366 ], shape = [9, 9], order = [2, 1])
367 cholref = reshape([real(TKG) :: 1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
368 , 1.0, 1.0, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
369 , 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM, DUM, DUM &
370 , 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM, DUM &
371 , 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM, DUM &
372 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM, DUM &
373 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM, DUM &
374 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, DUM &
375 , 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
376 ], shape = [9, 9], order = [2, 1])
377 call disp%skip()
378 call disp%show("cholmat")
379 call disp%show( cholmat , format = gform )
380 call disp%show("ndim = 9; roff = 0; coff = 0;")
381 ndim = 9; roff = 0; coff = 0;
382 if (present(bdim)) then
383 call disp%show("bdim")
384 call disp%show( bdim )
385 call disp%show("call setMatChol(cholmat, lowDia, info, iteration, ndim, roff, coff, bdim)")
386 call setMatChol(cholmat, lowDia, info, iteration, ndim, roff, coff, bdim)
387 else
388 call disp%show("call setMatChol(cholmat, lowDia, info, recursion, ndim, roff, coff)")
389 call setMatChol(cholmat, lowDia, info, recursion, ndim, roff, coff)
390 end if
391 call disp%show("cholmat")
392 call disp%show( cholmat , format = gform )
393 call disp%show("cholmat - cholref")
394 call disp%show( cholmat - cholref , format = gform )
395 call disp%show("info")
396 call disp%show( info )
397 call disp%skip()
398
399
400 cholmat = reshape([real(TKG) :: DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
401 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
402 , DUM, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
403 , DUM, DUM, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0 &
404 , DUM, DUM, DUM, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0 &
405 , DUM, DUM, DUM, DUM, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0 &
406 , DUM, DUM, DUM, DUM, DUM, 5.0, 5.0, 5.0, 5.0, 5.0 &
407 , DUM, DUM, DUM, DUM, DUM, DUM, 6.0, 6.0, 6.0, 6.0 &
408 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, 7.0, 7.0, 7.0 &
409 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 8.0, 8.0 &
410 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 9.0 &
411 ], shape = [11, 10], order = [2, 1])
412 cholref = reshape([real(TKG) :: DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
413 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM &
414 , DUM, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
415 , DUM, DUM, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
416 , DUM, DUM, DUM, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
417 , DUM, DUM, DUM, DUM, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 &
418 , DUM, DUM, DUM, DUM, DUM, 1.0, 1.0, 1.0, 1.0, 1.0 &
419 , DUM, DUM, DUM, DUM, DUM, DUM, 1.0, 1.0, 1.0, 1.0 &
420 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, 1.0, 1.0, 1.0 &
421 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 1.0, 1.0 &
422 , DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 1.0 &
423 ], shape = [11, 10], order = [2, 1])
424 call disp%skip()
425 call disp%show("cholmat")
426 call disp%show( cholmat , format = gform )
427 call disp%show("ndim = 9; roff = 2; coff = 1;")
428 ndim = 9; roff = 2; coff = 1;
429 if (present(bdim)) then
430 call disp%show("bdim")
431 call disp%show( bdim )
432 call disp%show("call setMatChol(cholmat, uppDia, info, iteration, ndim, roff, coff, bdim)")
433 call setMatChol(cholmat, uppDia, info, iteration, ndim, roff, coff, bdim)
434 else
435 call disp%show("call setMatChol(cholmat, uppDia, info, recursion, ndim, roff, coff)")
436 call setMatChol(cholmat, uppDia, info, recursion, ndim, roff, coff)
437 end if
438 call disp%show("cholmat")
439 call disp%show( cholmat , format = gform )
440 call disp%show("cholmat - cholref")
441 call disp%show( cholmat - cholref , format = gform )
442 call disp%show("info")
443 call disp%show( info )
444 call disp%skip()
445
446 end block
447
448 call disp%skip
449 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
450 call disp%show("!Cholesky factorization of complex positive-definite square matrix update.")
451 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
452 call disp%skip
453
454 block
455
456 complex(CKG), allocatable :: cholref(:,:), cholmat(:,:)
457
458 cholmat = reshape( [complex(CKG) :: CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM &
459 , CMPLX_DUMM, CMPLX_DUMM, (25.0, DUM), CMPLX_DUMM, CMPLX_DUMM &
460 , CMPLX_DUMM, CMPLX_DUMM, (-5.0, 5.0), (51.0, DUM), CMPLX_DUMM &
461 , CMPLX_DUMM, CMPLX_DUMM, (10.0, -5.0), ( 4.0, 6.0), (71.0, DUM) &
462 ], shape = [4, 5], order = [2, 1])
463 cholref = reshape( [complex(CKG) :: CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM &
464 , CMPLX_DUMM, CMPLX_DUMM, ( 5.0, 0.0), CMPLX_DUMM, CMPLX_DUMM &
465 , CMPLX_DUMM, CMPLX_DUMM, (-1.0, 1.0), (7.0, 0.0), CMPLX_DUMM &
466 , CMPLX_DUMM, CMPLX_DUMM, ( 2.0, -1.0), (1.0, 1.0), (8.0, 0.0) &
467 ], shape = [4, 5], order = [2, 1])
468 call disp%skip()
469 call disp%show("cholmat")
470 call disp%show( cholmat , format = cform )
471 call disp%show("ndim = 3; roff = 1; coff = 2;")
472 ndim = 3; roff = 1; coff = 2;
473 call disp%show("call setMatChol(cholmat, lowDia, info, iteration, ndim, roff, coff)")
474 call setMatChol(cholmat, lowDia, info, iteration, ndim, roff, coff)
475 call disp%show("cholmat")
476 call disp%show( cholmat , format = cform )
477 call disp%show("cholmat - cholref")
478 call disp%show( cholmat - cholref , format = cform )
479 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
480 if (info /= 0) error stop 'Cholesky factorization failed.'
481 call disp%skip()
482
483 cholmat = reshape( [complex(CKG) :: CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM &
484 , (9.0, DUM), ( 3.0,3.0), ( 3.0,-3.0), CMPLX_DUMM &
485 , CMPLX_DUMM, (18.0,DUM), ( 8.0,-6.0), CMPLX_DUMM &
486 , CMPLX_DUMM, CMPLX_DUMM, (43.0, DUM), CMPLX_DUMM &
487 ], shape = [4, 4], order = [2, 1])
488 cholref = reshape( [complex(CKG) :: CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM, CMPLX_DUMM &
489 , (3.0, 0.0), (1.0, 1.0), (1.0, -1.0), CMPLX_DUMM &
490 , CMPLX_DUMM, (4.0, 0.0), (2.0, -1.0), CMPLX_DUMM &
491 , CMPLX_DUMM, CMPLX_DUMM, (6.0, 0.0), CMPLX_DUMM &
492 ], shape = [4, 4], order = [2, 1])
493 call disp%skip()
494 call disp%show("cholmat")
495 call disp%show( cholmat , format = cform )
496 call disp%show("ndim = 3; roff = 1; coff = 0;")
497 ndim = 3; roff = 1; coff = 0;
498 if (present(bdim)) then
499 call disp%show("bdim")
500 call disp%show( bdim )
501 call disp%show("call setMatChol(cholmat, uppDia, info, iteration, ndim, roff, coff, bdim)")
502 call setMatChol(cholmat, uppDia, info, iteration, ndim, roff, coff, bdim)
503 else
504 call disp%show("call setMatChol(cholmat, uppDia, info, recursion, ndim, roff, coff)")
505 call setMatChol(cholmat, uppDia, info, recursion, ndim, roff, coff)
506 end if
507 call disp%show("cholmat")
508 call disp%show( cholmat , format = cform )
509 call disp%show("cholmat - cholref")
510 call disp%show( cholmat - cholref , format = cform )
511 call disp%show("if (info /= 0) error stop 'Cholesky factorization failed.'")
512 if (info /= 0) error stop 'Cholesky factorization failed.'
513 call disp%skip()
514
515 end block
516
517 end subroutine
518
519end program example
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Generate and return a generic or type/kind-specific IO format with the requested specifications that ...
Definition: pm_io.F90:18485
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter CKS
The single-precision complex kind in Fortran mode. On most platforms, this is a 32-bit real kind.
Definition: pm_kind.F90:570
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output
1
2mat
3+1.00000000, +0.00000000, +2.00000000
4+0.00000000, +4.00000000, +0.00000000
5+2.00000000, +0.00000000, +8.00000000
6choupp = 0
7call setMatChol(mat, uppDia, info, choupp, nothing)
8if (info /= 0) error stop
9choupp
10+1.00000000, +0.00000000, +2.00000000
11+0.00000000, +2.00000000, +0.00000000
12+0.00000000, +0.00000000, +2.00000000
13cholow = 0
14call setMatChol(mat, uppDia, info, cholow, transHerm)
15if (info /= 0) error stop
16cholow
17+1.00000000, +0.00000000, +0.00000000
18+0.00000000, +2.00000000, +0.00000000
19+2.00000000, +0.00000000, +2.00000000
20matmul(cholow, choupp) - mat ! must be all zero-valued elements.
21+0.00000000, +0.00000000, +0.00000000
22+0.00000000, +0.00000000, +0.00000000
23+0.00000000, +0.00000000, +0.00000000
24
25
26mat
27(+9.00000000, +0.00000000), (+3.00000000, +3.00000000), (+3.00000000, -3.00000000)
28(+3.00000000, -3.00000000), (+18.0000000, +0.00000000), (+8.00000000, -6.00000000)
29(+3.00000000, +3.00000000), (+8.00000000, +6.00000000), (+43.0000000, +0.00000000)
30choupp = 0
31call setMatChol(mat, uppDia, info, choupp, nothing)
32if (info /= 0) error stop
33choupp
34(+3.00000000, +0.00000000), (+1.00000000, +1.00000000), (+1.00000000, -1.00000000)
35(+0.00000000, +0.00000000), (+4.00000000, +0.00000000), (+2.00000000, -1.00000000)
36(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+6.00000000, +0.00000000)
37cholow = 0
38call setMatChol(mat, uppDia, info, cholow, transHerm)
39if (info /= 0) error stop
40cholow
41(+3.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
42(+1.00000000, -1.00000000), (+4.00000000, +0.00000000), (+0.00000000, +0.00000000)
43(+1.00000000, +1.00000000), (+2.00000000, +1.00000000), (+6.00000000, +0.00000000)
44matmul(cholow, choupp) - mat ! must be all zero-valued elements.
45(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
46(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
47(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
48
49
50mat
51(+25.0000000, +0.00000000), (-5.00000000, -5.00000000), (+10.0000000, +5.00000000)
52(-5.00000000, +5.00000000), (+51.0000000, +0.00000000), (+4.00000000, -6.00000000)
53(+10.0000000, -5.00000000), (+4.00000000, +6.00000000), (+71.0000000, +0.00000000)
54choupp = 0
55call setMatChol(mat, uppDia, info, choupp, nothing)
56if (info /= 0) error stop
57choupp
58(+5.00000000, +0.00000000), (-1.00000000, -1.00000000), (+2.00000000, +1.00000000)
59(+0.00000000, +0.00000000), (+7.00000000, +0.00000000), (+1.00000000, -1.00000000)
60(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+8.00000000, +0.00000000)
61cholow = 0
62call setMatChol(mat, uppDia, info, cholow, transHerm)
63if (info /= 0) error stop
64cholow
65(+5.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
66(-1.00000000, +1.00000000), (+7.00000000, +0.00000000), (+0.00000000, +0.00000000)
67(+2.00000000, -1.00000000), (+1.00000000, +1.00000000), (+8.00000000, +0.00000000)
68matmul(cholow, choupp) - mat ! must be all zero-valued elements.
69(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
70(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
71(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
72
73
74mat
75+1.00000000, +0.00000000, +2.00000000
76+0.00000000, +4.00000000, +0.00000000
77+2.00000000, +0.00000000, +8.00000000
78cholow = 0
79call setMatChol(mat, lowDia, info, cholow, nothing)
80if (info /= 0) error stop
81cholow
82+1.00000000, +0.00000000, +0.00000000
83+0.00000000, +2.00000000, +0.00000000
84+2.00000000, +0.00000000, +2.00000000
85choupp = 0
86choupp = 0
87call setMatChol(mat, lowDia, info, choupp, transHerm)
88if (info /= 0) error stop
89choupp
90+1.00000000, +0.00000000, +2.00000000
91+0.00000000, +2.00000000, +0.00000000
92+0.00000000, +0.00000000, +2.00000000
93matmul(cholow, choupp) - mat ! must be all zero-valued elements.
94+0.00000000, +0.00000000, +0.00000000
95+0.00000000, +0.00000000, +0.00000000
96+0.00000000, +0.00000000, +0.00000000
97
98
99mat
100(+9.00000000, +0.00000000), (+3.00000000, +3.00000000), (+3.00000000, -3.00000000)
101(+3.00000000, -3.00000000), (+18.0000000, +0.00000000), (+8.00000000, -6.00000000)
102(+3.00000000, +3.00000000), (+8.00000000, +6.00000000), (+43.0000000, +0.00000000)
103cholow = 0
104call setMatChol(mat, lowDia, info, cholow, nothing)
105if (info /= 0) error stop
106cholow
107(+3.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
108(+1.00000000, -1.00000000), (+4.00000000, +0.00000000), (+0.00000000, +0.00000000)
109(+1.00000000, +1.00000000), (+2.00000000, +1.00000000), (+6.00000000, +0.00000000)
110choupp = 0
111choupp = 0
112call setMatChol(mat, lowDia, info, choupp, transHerm)
113if (info /= 0) error stop
114choupp
115(+3.00000000, +0.00000000), (+1.00000000, +1.00000000), (+1.00000000, -1.00000000)
116(+0.00000000, +0.00000000), (+4.00000000, +0.00000000), (+2.00000000, -1.00000000)
117(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+6.00000000, +0.00000000)
118matmul(cholow, choupp) - mat ! must be all zero-valued elements.
119(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
120(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
121(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
122
123
124mat
125(+25.0000000, +0.00000000), (-5.00000000, -5.00000000), (+10.0000000, +5.00000000)
126(-5.00000000, +5.00000000), (+51.0000000, +0.00000000), (+4.00000000, -6.00000000)
127(+10.0000000, -5.00000000), (+4.00000000, +6.00000000), (+71.0000000, +0.00000000)
128cholow = 0
129call setMatChol(mat, lowDia, info, cholow, nothing)
130if (info /= 0) error stop
131cholow
132(+5.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
133(-1.00000000, +1.00000000), (+7.00000000, +0.00000000), (+0.00000000, +0.00000000)
134(+2.00000000, -1.00000000), (+1.00000000, +1.00000000), (+8.00000000, +0.00000000)
135choupp = 0
136choupp = 0
137call setMatChol(mat, lowDia, info, choupp, transHerm)
138if (info /= 0) error stop
139choupp
140(+5.00000000, +0.00000000), (-1.00000000, -1.00000000), (+2.00000000, +1.00000000)
141(+0.00000000, +0.00000000), (+7.00000000, +0.00000000), (+1.00000000, -1.00000000)
142(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+8.00000000, +0.00000000)
143matmul(cholow, choupp) - mat ! must be all zero-valued elements.
144(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
145(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
146(+0.00000000, +0.00000000), (+0.00000000, +0.00000000), (+0.00000000, +0.00000000)
147
148
149!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
150!In-place and out-of-place factorization within the same matrix
151!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152
153
154mat
155 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
156 +1.000000, +2.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
157 +1.000000, +2.000000, +3.000000, ***************, ***************, ***************, ***************, ***************, ***************
158 +1.000000, +2.000000, +3.000000, +4.000000, ***************, ***************, ***************, ***************, ***************
159 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, ***************, ***************, ***************, ***************
160 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, ***************, ***************, ***************
161 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, ***************, ***************
162 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, ***************
163 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +9.000000
164matchol(:,1:ndim) = mat
165matchol(:,1:ndim)
166 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
167 +1.000000, +2.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
168 +1.000000, +2.000000, +3.000000, ***************, ***************, ***************, ***************, ***************, ***************
169 +1.000000, +2.000000, +3.000000, +4.000000, ***************, ***************, ***************, ***************, ***************
170 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, ***************, ***************, ***************, ***************
171 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, ***************, ***************, ***************
172 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, ***************, ***************
173 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, ***************
174 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +9.000000
175call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)
176matchol(:,1:ndim)
177 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
178 +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
179 +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************, ***************
180 +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************
181 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************
182 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************
183 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************
184 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************
185 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
186matchol(:,1:ndim) - cholref
187 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
188 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
189 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
190 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
191 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
192 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
193 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
194 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
195 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
196if (info /= 0) error stop 'Cholesky factorization failed.'
197
198matchol(:,1:ndim) = mat
199matchol(:,1:ndim)
200 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
201 +1.000000, +2.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
202 +1.000000, +2.000000, +3.000000, ***************, ***************, ***************, ***************, ***************, ***************
203 +1.000000, +2.000000, +3.000000, +4.000000, ***************, ***************, ***************, ***************, ***************
204 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, ***************, ***************, ***************, ***************
205 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, ***************, ***************, ***************
206 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, ***************, ***************
207 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, ***************
208 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +9.000000
209call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)
210matchol(:,2:ndim+1)
211 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
212 +2.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
213 +2.000000, +3.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
214 +2.000000, +3.000000, +4.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
215 +2.000000, +3.000000, +4.000000, +5.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
216 +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +1.000000, +1.000000, +1.000000, +1.000000
217 +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +1.000000, +1.000000, +1.000000
218 +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +1.000000, +1.000000
219 +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +9.000000, +1.000000
220matchol(:,2:ndim+1) - transpose(cholref)
221 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
222***************, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
223***************, ***************, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
224***************, ***************, ***************, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
225***************, ***************, ***************, ***************, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
226***************, ***************, ***************, ***************, ***************, +0.000000, +0.000000, +0.000000, +0.000000
227***************, ***************, ***************, ***************, ***************, ***************, +0.000000, +0.000000, +0.000000
228***************, ***************, ***************, ***************, ***************, ***************, ***************, +0.000000, +0.000000
229***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, +0.000000
230if (info /= 0) error stop 'Cholesky factorization failed.'
231
232
233mat
234( +25.000000, +0.000000), (***************, ***************), (***************, ***************)
235( -5.000000, +5.000000), ( +51.000000, +0.000000), (***************, ***************)
236( +10.000000, -5.000000), ( +4.000000, +6.000000), ( +71.000000, +0.000000)
237matchol(:,1:ndim) = mat
238matchol(:,1:ndim)
239( +25.000000, +0.000000), (***************, ***************), (***************, ***************)
240( -5.000000, +5.000000), ( +51.000000, +0.000000), (***************, ***************)
241( +10.000000, -5.000000), ( +4.000000, +6.000000), ( +71.000000, +0.000000)
242call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,1:ndim), nothing)
243matchol(:,1:ndim)
244( +5.000000, +0.000000), (***************, ***************), (***************, ***************)
245( -1.000000, +1.000000), ( +7.000000, +0.000000), (***************, ***************)
246( +2.000000, -1.000000), ( +1.000000, +1.000000), ( +8.000000, +0.000000)
247cholref
248( +5.000000, +0.000000), (***************, ***************), (***************, ***************)
249( -1.000000, +1.000000), ( +7.000000, +0.000000), (***************, ***************)
250( +2.000000, -1.000000), ( +1.000000, +1.000000), ( +8.000000, +0.000000)
251matchol(:,1:ndim) - cholref
252( +0.000000, +0.000000), ( +0.000000, +0.000000), ( +0.000000, +0.000000)
253( +0.000000, +0.000000), ( +0.000000, +0.000000), ( +0.000000, +0.000000)
254( +0.000000, +0.000000), ( +0.000000, +0.000000), ( +0.000000, +0.000000)
255if (info /= 0) error stop 'Cholesky factorization failed.'
256
257matchol(:,1:ndim) = mat
258matchol(:,1:ndim)
259( +25.000000, +0.000000), (***************, ***************), (***************, ***************)
260( -5.000000, +5.000000), ( +51.000000, +0.000000), (***************, ***************)
261( +10.000000, -5.000000), ( +4.000000, +6.000000), ( +71.000000, +0.000000)
262call setMatChol(matchol(:,1:ndim), lowDia, info, matchol(:,2:ndim+1), transHerm)
263matchol(:,2:ndim+1)
264( +5.000000, +0.000000), ( -1.000000, -1.000000), ( +2.000000, +1.000000)
265( +51.000000, +0.000000), ( +7.000000, +0.000000), ( +1.000000, -1.000000)
266( +4.000000, +6.000000), ( +71.000000, +0.000000), ( +8.000000, +0.000000)
267transpose(conjg(cholref))
268( +5.000000, -0.000000), ( -1.000000, -1.000000), ( +2.000000, +1.000000)
269(***************, ***************), ( +7.000000, -0.000000), ( +1.000000, -1.000000)
270(***************, ***************), (***************, ***************), ( +8.000000, -0.000000)
271matchol(:,2:ndim+1) - transpose(conjg(cholref))
272( +0.000000, +0.000000), ( +0.000000, +0.000000), ( +0.000000, +0.000000)
273(***************, ***************), ( +0.000000, +0.000000), ( +0.000000, +0.000000)
274(***************, ***************), (***************, ***************), ( +0.000000, +0.000000)
275if (info /= 0) error stop 'Cholesky factorization failed.'
276
277
278!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
279!Cholesky factorization of real positive-definite square matrix update.
280!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281
282
283cholmat
284 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
285 +1.000000, +2.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
286 +1.000000, +2.000000, +3.000000, ***************, ***************, ***************, ***************, ***************, ***************
287 +1.000000, +2.000000, +3.000000, +4.000000, ***************, ***************, ***************, ***************, ***************
288 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, ***************, ***************, ***************, ***************
289 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, ***************, ***************, ***************
290 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, ***************, ***************
291 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, ***************
292 +1.000000, +2.000000, +3.000000, +4.000000, +5.000000, +6.000000, +7.000000, +8.000000, +9.000000
293ndim = 9; roff = 0; coff = 0;
294call setMatChol(cholmat, lowDia, info, recursion, ndim, roff, coff)
295cholmat
296 +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
297 +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************, ***************, ***************
298 +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************, ***************
299 +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************, ***************
300 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************, ***************
301 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************, ***************
302 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************, ***************
303 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, ***************
304 +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
305cholmat - cholref
306 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
307 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
308 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
309 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
310 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
311 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
312 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
313 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
314 +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000, +0.000000
315info
316+0
317
318
319cholmat
320***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
321***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
322***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
323***************, ***************, +2.000000, +2.000000, +2.000000, +2.000000, +2.000000, +2.000000, +2.000000, +2.000000
324***************, ***************, ***************, +3.000000, +3.000000, +3.000000, +3.000000, +3.000000, +3.000000, +3.000000
325***************, ***************, ***************, ***************, +4.000000, +4.000000, +4.000000, +4.000000, +4.000000, +4.000000
326***************, ***************, ***************, ***************, ***************, +5.000000, +5.000000, +5.000000, +5.000000, +5.000000
327***************, ***************, ***************, ***************, ***************, ***************, +6.000000, +6.000000, +6.000000, +6.000000
328***************, ***************, ***************, ***************, ***************, ***************, ***************, +7.000000, +7.000000, +7.000000
329***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, +8.000000, +8.000000
330***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, +9.000000
331ndim = 9; roff = 2; coff = 1;
332call setMatChol(cholmat, uppDia, info, recursion, ndim, roff, coff)
333cholmat
334***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
335***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************
336***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
337***************, ***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
338***************, ***************, ***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
339***************, ***************, ***************, ***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
340***************, ***************, ***************, ***************, ***************, +1.000000, +1.000000, +1.000000, +1.000000, +1.000000
341***************, ***************, ***************, ***************, ***************, ***************, +1.000000, +1.000000, +1.000000, +1.000000
342***************, ***************, ***************, ***************, ***************, ***************, ***************, +1.000000, +1.000000, +1.000000
343***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, +1.000000, +1.000000
344***************, ***************, ***************, ***************, ***************, ***************, ***************, ***************, *****