00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #ifndef TEMPLATE_LAPACK_POTRF_HEADER
00036 #define TEMPLATE_LAPACK_POTRF_HEADER
00037
00038 #include "template_lapack_potf2.h"
00039
00040 template<class Treal>
00041 int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *
00042 lda, integer *info)
00043 {
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 integer c__1 = 1;
00103 integer c_n1 = -1;
00104 Treal c_b13 = -1.;
00105 Treal c_b14 = 1.;
00106
00107
00108 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00109
00110 integer j;
00111 logical upper;
00112 integer jb, nb;
00113 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00114
00115
00116 a_dim1 = *lda;
00117 a_offset = 1 + a_dim1 * 1;
00118 a -= a_offset;
00119
00120
00121 *info = 0;
00122 upper = template_blas_lsame(uplo, "U");
00123 if (! upper && ! template_blas_lsame(uplo, "L")) {
00124 *info = -1;
00125 } else if (*n < 0) {
00126 *info = -2;
00127 } else if (*lda < maxMACRO(1,*n)) {
00128 *info = -4;
00129 }
00130 if (*info != 0) {
00131 i__1 = -(*info);
00132 template_blas_erbla("POTRF ", &i__1);
00133 return 0;
00134 }
00135
00136
00137
00138 if (*n == 0) {
00139 return 0;
00140 }
00141
00142
00143
00144 nb = template_lapack_ilaenv(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00145 ftnlen)1);
00146 if (nb <= 1 || nb >= *n) {
00147
00148
00149
00150 template_lapack_potf2(uplo, n, &a[a_offset], lda, info);
00151 } else {
00152
00153
00154
00155 if (upper) {
00156
00157
00158
00159 i__1 = *n;
00160 i__2 = nb;
00161 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00162
00163
00164
00165
00166
00167 i__3 = nb, i__4 = *n - j + 1;
00168 jb = minMACRO(i__3,i__4);
00169 i__3 = j - 1;
00170 template_blas_syrk("Upper", "Transpose", &jb, &i__3, &c_b13, &a_ref(1, j),
00171 lda, &c_b14, &a_ref(j, j), lda)
00172 ;
00173 template_lapack_potf2("Upper", &jb, &a_ref(j, j), lda, info);
00174 if (*info != 0) {
00175 goto L30;
00176 }
00177 if (j + jb <= *n) {
00178
00179
00180
00181 i__3 = *n - j - jb + 1;
00182 i__4 = j - 1;
00183 template_blas_gemm("Transpose", "No transpose", &jb, &i__3, &i__4, &
00184 c_b13, &a_ref(1, j), lda, &a_ref(1, j + jb), lda,
00185 &c_b14, &a_ref(j, j + jb), lda);
00186 i__3 = *n - j - jb + 1;
00187 template_blas_trsm("Left", "Upper", "Transpose", "Non-unit", &jb, &
00188 i__3, &c_b14, &a_ref(j, j), lda, &a_ref(j, j + jb)
00189 , lda)
00190 ;
00191 }
00192
00193 }
00194
00195 } else {
00196
00197
00198
00199 i__2 = *n;
00200 i__1 = nb;
00201 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00202
00203
00204
00205
00206
00207 i__3 = nb, i__4 = *n - j + 1;
00208 jb = minMACRO(i__3,i__4);
00209 i__3 = j - 1;
00210 template_blas_syrk("Lower", "No transpose", &jb, &i__3, &c_b13, &a_ref(j,
00211 1), lda, &c_b14, &a_ref(j, j), lda);
00212 template_lapack_potf2("Lower", &jb, &a_ref(j, j), lda, info);
00213 if (*info != 0) {
00214 goto L30;
00215 }
00216 if (j + jb <= *n) {
00217
00218
00219
00220 i__3 = *n - j - jb + 1;
00221 i__4 = j - 1;
00222 template_blas_gemm("No transpose", "Transpose", &i__3, &jb, &i__4, &
00223 c_b13, &a_ref(j + jb, 1), lda, &a_ref(j, 1), lda,
00224 &c_b14, &a_ref(j + jb, j), lda);
00225 i__3 = *n - j - jb + 1;
00226 template_blas_trsm("Right", "Lower", "Transpose", "Non-unit", &i__3, &
00227 jb, &c_b14, &a_ref(j, j), lda, &a_ref(j + jb, j),
00228 lda);
00229 }
00230
00231 }
00232 }
00233 }
00234 goto L40;
00235
00236 L30:
00237 *info = *info + j - 1;
00238
00239 L40:
00240 return 0;
00241
00242
00243
00244 }
00245
00246 #undef a_ref
00247
00248
00249 #endif