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
00036
00037 #ifndef TEMPLATE_LAPACK_POTRF_HEADER
00038 #define TEMPLATE_LAPACK_POTRF_HEADER
00039
00040 #include "template_lapack_potf2.h"
00041
00042 template<class Treal>
00043 int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *
00044 lda, integer *info)
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
00103
00104 integer c__1 = 1;
00105 integer c_n1 = -1;
00106 Treal c_b13 = -1.;
00107 Treal c_b14 = 1.;
00108
00109
00110 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00111
00112 integer j;
00113 logical upper;
00114 integer jb, nb;
00115 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00116
00117
00118 a_dim1 = *lda;
00119 a_offset = 1 + a_dim1 * 1;
00120 a -= a_offset;
00121
00122
00123 *info = 0;
00124 upper = template_blas_lsame(uplo, "U");
00125 if (! upper && ! template_blas_lsame(uplo, "L")) {
00126 *info = -1;
00127 } else if (*n < 0) {
00128 *info = -2;
00129 } else if (*lda < maxMACRO(1,*n)) {
00130 *info = -4;
00131 }
00132 if (*info != 0) {
00133 i__1 = -(*info);
00134 template_blas_erbla("POTRF ", &i__1);
00135 return 0;
00136 }
00137
00138
00139
00140 if (*n == 0) {
00141 return 0;
00142 }
00143
00144
00145
00146 nb = template_lapack_ilaenv(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00147 ftnlen)1);
00148 if (nb <= 1 || nb >= *n) {
00149
00150
00151
00152 template_lapack_potf2(uplo, n, &a[a_offset], lda, info);
00153 } else {
00154
00155
00156
00157 if (upper) {
00158
00159
00160
00161 i__1 = *n;
00162 i__2 = nb;
00163 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00164
00165
00166
00167
00168
00169 i__3 = nb, i__4 = *n - j + 1;
00170 jb = minMACRO(i__3,i__4);
00171 i__3 = j - 1;
00172 template_blas_syrk("Upper", "Transpose", &jb, &i__3, &c_b13, &a_ref(1, j),
00173 lda, &c_b14, &a_ref(j, j), lda)
00174 ;
00175 template_lapack_potf2("Upper", &jb, &a_ref(j, j), lda, info);
00176 if (*info != 0) {
00177 goto L30;
00178 }
00179 if (j + jb <= *n) {
00180
00181
00182
00183 i__3 = *n - j - jb + 1;
00184 i__4 = j - 1;
00185 template_blas_gemm("Transpose", "No transpose", &jb, &i__3, &i__4, &
00186 c_b13, &a_ref(1, j), lda, &a_ref(1, j + jb), lda,
00187 &c_b14, &a_ref(j, j + jb), lda);
00188 i__3 = *n - j - jb + 1;
00189 template_blas_trsm("Left", "Upper", "Transpose", "Non-unit", &jb, &
00190 i__3, &c_b14, &a_ref(j, j), lda, &a_ref(j, j + jb)
00191 , lda)
00192 ;
00193 }
00194
00195 }
00196
00197 } else {
00198
00199
00200
00201 i__2 = *n;
00202 i__1 = nb;
00203 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00204
00205
00206
00207
00208
00209 i__3 = nb, i__4 = *n - j + 1;
00210 jb = minMACRO(i__3,i__4);
00211 i__3 = j - 1;
00212 template_blas_syrk("Lower", "No transpose", &jb, &i__3, &c_b13, &a_ref(j,
00213 1), lda, &c_b14, &a_ref(j, j), lda);
00214 template_lapack_potf2("Lower", &jb, &a_ref(j, j), lda, info);
00215 if (*info != 0) {
00216 goto L30;
00217 }
00218 if (j + jb <= *n) {
00219
00220
00221
00222 i__3 = *n - j - jb + 1;
00223 i__4 = j - 1;
00224 template_blas_gemm("No transpose", "Transpose", &i__3, &jb, &i__4, &
00225 c_b13, &a_ref(j + jb, 1), lda, &a_ref(j, 1), lda,
00226 &c_b14, &a_ref(j + jb, j), lda);
00227 i__3 = *n - j - jb + 1;
00228 template_blas_trsm("Right", "Lower", "Transpose", "Non-unit", &i__3, &
00229 jb, &c_b14, &a_ref(j, j), lda, &a_ref(j + jb, j),
00230 lda);
00231 }
00232
00233 }
00234 }
00235 }
00236 goto L40;
00237
00238 L30:
00239 *info = *info + j - 1;
00240
00241 L40:
00242 return 0;
00243
00244
00245
00246 }
00247
00248 #undef a_ref
00249
00250
00251 #endif