1 /* -- translated by f2c (version 20191129). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11 */ 12 13 #include "f2c.h" 14 15 /* Table of constant values */ 16 17 static integer c__1 = 1; 18 static integer c_n1 = -1; 19 static doublereal c_b13 = -1.; 20 static doublereal c_b14 = 1.; 21 22 /* > \brief \b DPOTRF 23 24 =========== DOCUMENTATION =========== 25 26 Online html documentation available at 27 http://www.netlib.org/lapack/explore-html/ 28 29 > \htmlonly 30 > Download DPOTRF + dependencies 31 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpotrf. 32 f"> 33 > [TGZ]</a> 34 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpotrf. 35 f"> 36 > [ZIP]</a> 37 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpotrf. 38 f"> 39 > [TXT]</a> 40 > \endhtmlonly 41 42 Definition: 43 =========== 44 45 SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO ) 46 47 CHARACTER UPLO 48 INTEGER INFO, LDA, N 49 DOUBLE PRECISION A( LDA, * ) 50 51 52 > \par Purpose: 53 ============= 54 > 55 > \verbatim 56 > 57 > DPOTRF computes the Cholesky factorization of a real symmetric 58 > positive definite matrix A. 59 > 60 > The factorization has the form 61 > A = U**T * U, if UPLO = 'U', or 62 > A = L * L**T, if UPLO = 'L', 63 > where U is an upper triangular matrix and L is lower triangular. 64 > 65 > This is the block version of the algorithm, calling Level 3 BLAS. 66 > \endverbatim 67 68 Arguments: 69 ========== 70 71 > \param[in] UPLO 72 > \verbatim 73 > UPLO is CHARACTER*1 74 > = 'U': Upper triangle of A is stored; 75 > = 'L': Lower triangle of A is stored. 76 > \endverbatim 77 > 78 > \param[in] N 79 > \verbatim 80 > N is INTEGER 81 > The order of the matrix A. N >= 0. 82 > \endverbatim 83 > 84 > \param[in,out] A 85 > \verbatim 86 > A is DOUBLE PRECISION array, dimension (LDA,N) 87 > On entry, the symmetric matrix A. If UPLO = 'U', the leading 88 > N-by-N upper triangular part of A contains the upper 89 > triangular part of the matrix A, and the strictly lower 90 > triangular part of A is not referenced. If UPLO = 'L', the 91 > leading N-by-N lower triangular part of A contains the lower 92 > triangular part of the matrix A, and the strictly upper 93 > triangular part of A is not referenced. 94 > 95 > On exit, if INFO = 0, the factor U or L from the Cholesky 96 > factorization A = U**T*U or A = L*L**T. 97 > \endverbatim 98 > 99 > \param[in] LDA 100 > \verbatim 101 > LDA is INTEGER 102 > The leading dimension of the array A. LDA >= max(1,N). 103 > \endverbatim 104 > 105 > \param[out] INFO 106 > \verbatim 107 > INFO is INTEGER 108 > = 0: successful exit 109 > < 0: if INFO = -i, the i-th argument had an illegal value 110 > > 0: if INFO = i, the leading minor of order i is not 111 > positive definite, and the factorization could not be 112 > completed. 113 > \endverbatim 114 115 Authors: 116 ======== 117 118 > \author Univ. of Tennessee 119 > \author Univ. of California Berkeley 120 > \author Univ. of Colorado Denver 121 > \author NAG Ltd. 122 123 > \date November 2011 124 125 > \ingroup doublePOcomputational 126 127 ===================================================================== igraphdpotrf_(char * uplo,integer * n,doublereal * a,integer * lda,integer * info)128 Subroutine */ int igraphdpotrf_(char *uplo, integer *n, doublereal *a, integer * 129 lda, integer *info) 130 { 131 /* System generated locals */ 132 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 133 134 /* Local variables */ 135 integer j, jb, nb; 136 extern /* Subroutine */ int igraphdgemm_(char *, char *, integer *, integer *, 137 integer *, doublereal *, doublereal *, integer *, doublereal *, 138 integer *, doublereal *, doublereal *, integer *); 139 extern logical igraphlsame_(char *, char *); 140 extern /* Subroutine */ int igraphdtrsm_(char *, char *, char *, char *, 141 integer *, integer *, doublereal *, doublereal *, integer *, 142 doublereal *, integer *); 143 logical upper; 144 extern /* Subroutine */ int igraphdsyrk_(char *, char *, integer *, integer *, 145 doublereal *, doublereal *, integer *, doublereal *, doublereal *, 146 integer *), igraphdpotf2_(char *, integer *, 147 doublereal *, integer *, integer *), igraphxerbla_(char *, 148 integer *, ftnlen); 149 extern integer igraphilaenv_(integer *, char *, char *, integer *, integer *, 150 integer *, integer *, ftnlen, ftnlen); 151 152 153 /* -- LAPACK computational routine (version 3.4.0) -- 154 -- LAPACK is a software package provided by Univ. of Tennessee, -- 155 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 156 November 2011 157 158 159 ===================================================================== 160 161 162 Test the input parameters. 163 164 Parameter adjustments */ 165 a_dim1 = *lda; 166 a_offset = 1 + a_dim1; 167 a -= a_offset; 168 169 /* Function Body */ 170 *info = 0; 171 upper = igraphlsame_(uplo, "U"); 172 if (! upper && ! igraphlsame_(uplo, "L")) { 173 *info = -1; 174 } else if (*n < 0) { 175 *info = -2; 176 } else if (*lda < max(1,*n)) { 177 *info = -4; 178 } 179 if (*info != 0) { 180 i__1 = -(*info); 181 igraphxerbla_("DPOTRF", &i__1, (ftnlen)6); 182 return 0; 183 } 184 185 /* Quick return if possible */ 186 187 if (*n == 0) { 188 return 0; 189 } 190 191 /* Determine the block size for this environment. */ 192 193 nb = igraphilaenv_(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( 194 ftnlen)1); 195 if (nb <= 1 || nb >= *n) { 196 197 /* Use unblocked code. */ 198 199 igraphdpotf2_(uplo, n, &a[a_offset], lda, info); 200 } else { 201 202 /* Use blocked code. */ 203 204 if (upper) { 205 206 /* Compute the Cholesky factorization A = U**T*U. */ 207 208 i__1 = *n; 209 i__2 = nb; 210 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { 211 212 /* Update and factorize the current diagonal block and test 213 for non-positive-definiteness. 214 215 Computing MIN */ 216 i__3 = nb, i__4 = *n - j + 1; 217 jb = min(i__3,i__4); 218 i__3 = j - 1; 219 igraphdsyrk_("Upper", "Transpose", &jb, &i__3, &c_b13, &a[j * 220 a_dim1 + 1], lda, &c_b14, &a[j + j * a_dim1], lda); 221 igraphdpotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info); 222 if (*info != 0) { 223 goto L30; 224 } 225 if (j + jb <= *n) { 226 227 /* Compute the current block row. */ 228 229 i__3 = *n - j - jb + 1; 230 i__4 = j - 1; 231 igraphdgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, & 232 c_b13, &a[j * a_dim1 + 1], lda, &a[(j + jb) * 233 a_dim1 + 1], lda, &c_b14, &a[j + (j + jb) * 234 a_dim1], lda); 235 i__3 = *n - j - jb + 1; 236 igraphdtrsm_("Left", "Upper", "Transpose", "Non-unit", &jb, & 237 i__3, &c_b14, &a[j + j * a_dim1], lda, &a[j + (j 238 + jb) * a_dim1], lda); 239 } 240 /* L10: */ 241 } 242 243 } else { 244 245 /* Compute the Cholesky factorization A = L*L**T. */ 246 247 i__2 = *n; 248 i__1 = nb; 249 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { 250 251 /* Update and factorize the current diagonal block and test 252 for non-positive-definiteness. 253 254 Computing MIN */ 255 i__3 = nb, i__4 = *n - j + 1; 256 jb = min(i__3,i__4); 257 i__3 = j - 1; 258 igraphdsyrk_("Lower", "No transpose", &jb, &i__3, &c_b13, &a[j + 259 a_dim1], lda, &c_b14, &a[j + j * a_dim1], lda); 260 igraphdpotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info); 261 if (*info != 0) { 262 goto L30; 263 } 264 if (j + jb <= *n) { 265 266 /* Compute the current block column. */ 267 268 i__3 = *n - j - jb + 1; 269 i__4 = j - 1; 270 igraphdgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, & 271 c_b13, &a[j + jb + a_dim1], lda, &a[j + a_dim1], 272 lda, &c_b14, &a[j + jb + j * a_dim1], lda); 273 i__3 = *n - j - jb + 1; 274 igraphdtrsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, & 275 jb, &c_b14, &a[j + j * a_dim1], lda, &a[j + jb + 276 j * a_dim1], lda); 277 } 278 /* L20: */ 279 } 280 } 281 } 282 goto L40; 283 284 L30: 285 *info = *info + j - 1; 286 287 L40: 288 return 0; 289 290 /* End of DPOTRF */ 291 292 } /* igraphdpotrf_ */ 293 294