ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_LARFB_HEADER 00036 #define TEMPLATE_LAPACK_LARFB_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_larfb(const char *side, const char *trans, const char *direct, const char * 00041 storev, const integer *m, const integer *n, const integer *k, const Treal *v, const integer * 00042 ldv, const Treal *t, const integer *ldt, Treal *c__, const integer *ldc, 00043 Treal *work, const integer *ldwork) 00044 { 00045 /* -- LAPACK auxiliary routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 February 29, 1992 00049 00050 00051 Purpose 00052 ======= 00053 00054 DLARFB applies a real block reflector H or its transpose H' to a 00055 real m by n matrix C, from either the left or the right. 00056 00057 Arguments 00058 ========= 00059 00060 SIDE (input) CHARACTER*1 00061 = 'L': apply H or H' from the Left 00062 = 'R': apply H or H' from the Right 00063 00064 TRANS (input) CHARACTER*1 00065 = 'N': apply H (No transpose) 00066 = 'T': apply H' (Transpose) 00067 00068 DIRECT (input) CHARACTER*1 00069 Indicates how H is formed from a product of elementary 00070 reflectors 00071 = 'F': H = H(1) H(2) . . . H(k) (Forward) 00072 = 'B': H = H(k) . . . H(2) H(1) (Backward) 00073 00074 STOREV (input) CHARACTER*1 00075 Indicates how the vectors which define the elementary 00076 reflectors are stored: 00077 = 'C': Columnwise 00078 = 'R': Rowwise 00079 00080 M (input) INTEGER 00081 The number of rows of the matrix C. 00082 00083 N (input) INTEGER 00084 The number of columns of the matrix C. 00085 00086 K (input) INTEGER 00087 The order of the matrix T (= the number of elementary 00088 reflectors whose product defines the block reflector). 00089 00090 V (input) DOUBLE PRECISION array, dimension 00091 (LDV,K) if STOREV = 'C' 00092 (LDV,M) if STOREV = 'R' and SIDE = 'L' 00093 (LDV,N) if STOREV = 'R' and SIDE = 'R' 00094 The matrix V. See further details. 00095 00096 LDV (input) INTEGER 00097 The leading dimension of the array V. 00098 If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); 00099 if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); 00100 if STOREV = 'R', LDV >= K. 00101 00102 T (input) DOUBLE PRECISION array, dimension (LDT,K) 00103 The triangular k by k matrix T in the representation of the 00104 block reflector. 00105 00106 LDT (input) INTEGER 00107 The leading dimension of the array T. LDT >= K. 00108 00109 C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00110 On entry, the m by n matrix C. 00111 On exit, C is overwritten by H*C or H'*C or C*H or C*H'. 00112 00113 LDC (input) INTEGER 00114 The leading dimension of the array C. LDA >= max(1,M). 00115 00116 WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) 00117 00118 LDWORK (input) INTEGER 00119 The leading dimension of the array WORK. 00120 If SIDE = 'L', LDWORK >= max(1,N); 00121 if SIDE = 'R', LDWORK >= max(1,M). 00122 00123 ===================================================================== 00124 00125 00126 Quick return if possible 00127 00128 Parameter adjustments */ 00129 /* Table of constant values */ 00130 integer c__1 = 1; 00131 Treal c_b14 = 1.; 00132 Treal c_b25 = -1.; 00133 00134 /* System generated locals */ 00135 integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1, 00136 work_offset, i__1, i__2; 00137 /* Local variables */ 00138 integer i__, j; 00139 char transt[1]; 00140 #define work_ref(a_1,a_2) work[(a_2)*work_dim1 + a_1] 00141 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1] 00142 #define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1] 00143 00144 00145 v_dim1 = *ldv; 00146 v_offset = 1 + v_dim1 * 1; 00147 v -= v_offset; 00148 t_dim1 = *ldt; 00149 t_offset = 1 + t_dim1 * 1; 00150 t -= t_offset; 00151 c_dim1 = *ldc; 00152 c_offset = 1 + c_dim1 * 1; 00153 c__ -= c_offset; 00154 work_dim1 = *ldwork; 00155 work_offset = 1 + work_dim1 * 1; 00156 work -= work_offset; 00157 00158 /* Function Body */ 00159 if (*m <= 0 || *n <= 0) { 00160 return 0; 00161 } 00162 00163 if (template_blas_lsame(trans, "N")) { 00164 *(unsigned char *)transt = 'T'; 00165 } else { 00166 *(unsigned char *)transt = 'N'; 00167 } 00168 00169 if (template_blas_lsame(storev, "C")) { 00170 00171 if (template_blas_lsame(direct, "F")) { 00172 00173 /* Let V = ( V1 ) (first K rows) 00174 ( V2 ) 00175 where V1 is unit lower triangular. */ 00176 00177 if (template_blas_lsame(side, "L")) { 00178 00179 /* Form H * C or H' * C where C = ( C1 ) 00180 ( C2 ) 00181 00182 W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 00183 00184 W := C1' */ 00185 00186 i__1 = *k; 00187 for (j = 1; j <= i__1; ++j) { 00188 template_blas_copy(n, &c___ref(j, 1), ldc, &work_ref(1, j), &c__1); 00189 /* L10: */ 00190 } 00191 00192 /* W := W * V1 */ 00193 00194 template_blas_trmm("Right", "Lower", "No transpose", "Unit", n, k, &c_b14, 00195 &v[v_offset], ldv, &work[work_offset], ldwork); 00196 if (*m > *k) { 00197 00198 /* W := W + C2'*V2 */ 00199 00200 i__1 = *m - *k; 00201 template_blas_gemm("Transpose", "No transpose", n, k, &i__1, &c_b14, & 00202 c___ref(*k + 1, 1), ldc, &v_ref(*k + 1, 1), ldv, & 00203 c_b14, &work[work_offset], ldwork); 00204 } 00205 00206 /* W := W * T' or W * T */ 00207 00208 template_blas_trmm("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[ 00209 t_offset], ldt, &work[work_offset], ldwork); 00210 00211 /* C := C - V * W' */ 00212 00213 if (*m > *k) { 00214 00215 /* C2 := C2 - V2 * W' */ 00216 00217 i__1 = *m - *k; 00218 template_blas_gemm("No transpose", "Transpose", &i__1, n, k, &c_b25, & 00219 v_ref(*k + 1, 1), ldv, &work[work_offset], ldwork, 00220 &c_b14, &c___ref(*k + 1, 1), ldc); 00221 } 00222 00223 /* W := W * V1' */ 00224 00225 template_blas_trmm("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, & 00226 v[v_offset], ldv, &work[work_offset], ldwork); 00227 00228 /* C1 := C1 - W' */ 00229 00230 i__1 = *k; 00231 for (j = 1; j <= i__1; ++j) { 00232 i__2 = *n; 00233 for (i__ = 1; i__ <= i__2; ++i__) { 00234 c___ref(j, i__) = c___ref(j, i__) - work_ref(i__, j); 00235 /* L20: */ 00236 } 00237 /* L30: */ 00238 } 00239 00240 } else if (template_blas_lsame(side, "R")) { 00241 00242 /* Form C * H or C * H' where C = ( C1 C2 ) 00243 00244 W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00245 00246 W := C1 */ 00247 00248 i__1 = *k; 00249 for (j = 1; j <= i__1; ++j) { 00250 template_blas_copy(m, &c___ref(1, j), &c__1, &work_ref(1, j), &c__1); 00251 /* L40: */ 00252 } 00253 00254 /* W := W * V1 */ 00255 00256 template_blas_trmm("Right", "Lower", "No transpose", "Unit", m, k, &c_b14, 00257 &v[v_offset], ldv, &work[work_offset], ldwork); 00258 if (*n > *k) { 00259 00260 /* W := W + C2 * V2 */ 00261 00262 i__1 = *n - *k; 00263 template_blas_gemm("No transpose", "No transpose", m, k, &i__1, & 00264 c_b14, &c___ref(1, *k + 1), ldc, &v_ref(*k + 1, 1) 00265 , ldv, &c_b14, &work[work_offset], ldwork); 00266 } 00267 00268 /* W := W * T or W * T' */ 00269 00270 template_blas_trmm("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[ 00271 t_offset], ldt, &work[work_offset], ldwork); 00272 00273 /* C := C - W * V' */ 00274 00275 if (*n > *k) { 00276 00277 /* C2 := C2 - W * V2' */ 00278 00279 i__1 = *n - *k; 00280 template_blas_gemm("No transpose", "Transpose", m, &i__1, k, &c_b25, & 00281 work[work_offset], ldwork, &v_ref(*k + 1, 1), ldv, 00282 &c_b14, &c___ref(1, *k + 1), ldc); 00283 } 00284 00285 /* W := W * V1' */ 00286 00287 template_blas_trmm("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, & 00288 v[v_offset], ldv, &work[work_offset], ldwork); 00289 00290 /* C1 := C1 - W */ 00291 00292 i__1 = *k; 00293 for (j = 1; j <= i__1; ++j) { 00294 i__2 = *m; 00295 for (i__ = 1; i__ <= i__2; ++i__) { 00296 c___ref(i__, j) = c___ref(i__, j) - work_ref(i__, j); 00297 /* L50: */ 00298 } 00299 /* L60: */ 00300 } 00301 } 00302 00303 } else { 00304 00305 /* Let V = ( V1 ) 00306 ( V2 ) (last K rows) 00307 where V2 is unit upper triangular. */ 00308 00309 if (template_blas_lsame(side, "L")) { 00310 00311 /* Form H * C or H' * C where C = ( C1 ) 00312 ( C2 ) 00313 00314 W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 00315 00316 W := C2' */ 00317 00318 i__1 = *k; 00319 for (j = 1; j <= i__1; ++j) { 00320 template_blas_copy(n, &c___ref(*m - *k + j, 1), ldc, &work_ref(1, j), 00321 &c__1); 00322 /* L70: */ 00323 } 00324 00325 /* W := W * V2 */ 00326 00327 template_blas_trmm("Right", "Upper", "No transpose", "Unit", n, k, &c_b14, 00328 &v_ref(*m - *k + 1, 1), ldv, &work[work_offset], 00329 ldwork); 00330 if (*m > *k) { 00331 00332 /* W := W + C1'*V1 */ 00333 00334 i__1 = *m - *k; 00335 template_blas_gemm("Transpose", "No transpose", n, k, &i__1, &c_b14, & 00336 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, & 00337 work[work_offset], ldwork); 00338 } 00339 00340 /* W := W * T' or W * T */ 00341 00342 template_blas_trmm("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[ 00343 t_offset], ldt, &work[work_offset], ldwork); 00344 00345 /* C := C - V * W' */ 00346 00347 if (*m > *k) { 00348 00349 /* C1 := C1 - V1 * W' */ 00350 00351 i__1 = *m - *k; 00352 template_blas_gemm("No transpose", "Transpose", &i__1, n, k, &c_b25, & 00353 v[v_offset], ldv, &work[work_offset], ldwork, & 00354 c_b14, &c__[c_offset], ldc) 00355 ; 00356 } 00357 00358 /* W := W * V2' */ 00359 00360 template_blas_trmm("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, & 00361 v_ref(*m - *k + 1, 1), ldv, &work[work_offset], 00362 ldwork); 00363 00364 /* C2 := C2 - W' */ 00365 00366 i__1 = *k; 00367 for (j = 1; j <= i__1; ++j) { 00368 i__2 = *n; 00369 for (i__ = 1; i__ <= i__2; ++i__) { 00370 c___ref(*m - *k + j, i__) = c___ref(*m - *k + j, i__) 00371 - work_ref(i__, j); 00372 /* L80: */ 00373 } 00374 /* L90: */ 00375 } 00376 00377 } else if (template_blas_lsame(side, "R")) { 00378 00379 /* Form C * H or C * H' where C = ( C1 C2 ) 00380 00381 W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00382 00383 W := C2 */ 00384 00385 i__1 = *k; 00386 for (j = 1; j <= i__1; ++j) { 00387 template_blas_copy(m, &c___ref(1, *n - *k + j), &c__1, &work_ref(1, j) 00388 , &c__1); 00389 /* L100: */ 00390 } 00391 00392 /* W := W * V2 */ 00393 00394 template_blas_trmm("Right", "Upper", "No transpose", "Unit", m, k, &c_b14, 00395 &v_ref(*n - *k + 1, 1), ldv, &work[work_offset], 00396 ldwork); 00397 if (*n > *k) { 00398 00399 /* W := W + C1 * V1 */ 00400 00401 i__1 = *n - *k; 00402 template_blas_gemm("No transpose", "No transpose", m, k, &i__1, & 00403 c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, & 00404 c_b14, &work[work_offset], ldwork); 00405 } 00406 00407 /* W := W * T or W * T' */ 00408 00409 template_blas_trmm("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[ 00410 t_offset], ldt, &work[work_offset], ldwork); 00411 00412 /* C := C - W * V' */ 00413 00414 if (*n > *k) { 00415 00416 /* C1 := C1 - W * V1' */ 00417 00418 i__1 = *n - *k; 00419 template_blas_gemm("No transpose", "Transpose", m, &i__1, k, &c_b25, & 00420 work[work_offset], ldwork, &v[v_offset], ldv, & 00421 c_b14, &c__[c_offset], ldc) 00422 ; 00423 } 00424 00425 /* W := W * V2' */ 00426 00427 template_blas_trmm("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, & 00428 v_ref(*n - *k + 1, 1), ldv, &work[work_offset], 00429 ldwork); 00430 00431 /* C2 := C2 - W */ 00432 00433 i__1 = *k; 00434 for (j = 1; j <= i__1; ++j) { 00435 i__2 = *m; 00436 for (i__ = 1; i__ <= i__2; ++i__) { 00437 c___ref(i__, *n - *k + j) = c___ref(i__, *n - *k + j) 00438 - work_ref(i__, j); 00439 /* L110: */ 00440 } 00441 /* L120: */ 00442 } 00443 } 00444 } 00445 00446 } else if (template_blas_lsame(storev, "R")) { 00447 00448 if (template_blas_lsame(direct, "F")) { 00449 00450 /* Let V = ( V1 V2 ) (V1: first K columns) 00451 where V1 is unit upper triangular. */ 00452 00453 if (template_blas_lsame(side, "L")) { 00454 00455 /* Form H * C or H' * C where C = ( C1 ) 00456 ( C2 ) 00457 00458 W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 00459 00460 W := C1' */ 00461 00462 i__1 = *k; 00463 for (j = 1; j <= i__1; ++j) { 00464 template_blas_copy(n, &c___ref(j, 1), ldc, &work_ref(1, j), &c__1); 00465 /* L130: */ 00466 } 00467 00468 /* W := W * V1' */ 00469 00470 template_blas_trmm("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, & 00471 v[v_offset], ldv, &work[work_offset], ldwork); 00472 if (*m > *k) { 00473 00474 /* W := W + C2'*V2' */ 00475 00476 i__1 = *m - *k; 00477 template_blas_gemm("Transpose", "Transpose", n, k, &i__1, &c_b14, & 00478 c___ref(*k + 1, 1), ldc, &v_ref(1, *k + 1), ldv, & 00479 c_b14, &work[work_offset], ldwork); 00480 } 00481 00482 /* W := W * T' or W * T */ 00483 00484 template_blas_trmm("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[ 00485 t_offset], ldt, &work[work_offset], ldwork); 00486 00487 /* C := C - V' * W' */ 00488 00489 if (*m > *k) { 00490 00491 /* C2 := C2 - V2' * W' */ 00492 00493 i__1 = *m - *k; 00494 template_blas_gemm("Transpose", "Transpose", &i__1, n, k, &c_b25, & 00495 v_ref(1, *k + 1), ldv, &work[work_offset], ldwork, 00496 &c_b14, &c___ref(*k + 1, 1), ldc); 00497 } 00498 00499 /* W := W * V1 */ 00500 00501 template_blas_trmm("Right", "Upper", "No transpose", "Unit", n, k, &c_b14, 00502 &v[v_offset], ldv, &work[work_offset], ldwork); 00503 00504 /* C1 := C1 - W' */ 00505 00506 i__1 = *k; 00507 for (j = 1; j <= i__1; ++j) { 00508 i__2 = *n; 00509 for (i__ = 1; i__ <= i__2; ++i__) { 00510 c___ref(j, i__) = c___ref(j, i__) - work_ref(i__, j); 00511 /* L140: */ 00512 } 00513 /* L150: */ 00514 } 00515 00516 } else if (template_blas_lsame(side, "R")) { 00517 00518 /* Form C * H or C * H' where C = ( C1 C2 ) 00519 00520 W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 00521 00522 W := C1 */ 00523 00524 i__1 = *k; 00525 for (j = 1; j <= i__1; ++j) { 00526 template_blas_copy(m, &c___ref(1, j), &c__1, &work_ref(1, j), &c__1); 00527 /* L160: */ 00528 } 00529 00530 /* W := W * V1' */ 00531 00532 template_blas_trmm("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, & 00533 v[v_offset], ldv, &work[work_offset], ldwork); 00534 if (*n > *k) { 00535 00536 /* W := W + C2 * V2' */ 00537 00538 i__1 = *n - *k; 00539 template_blas_gemm("No transpose", "Transpose", m, k, &i__1, &c_b14, & 00540 c___ref(1, *k + 1), ldc, &v_ref(1, *k + 1), ldv, & 00541 c_b14, &work[work_offset], ldwork); 00542 } 00543 00544 /* W := W * T or W * T' */ 00545 00546 template_blas_trmm("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[ 00547 t_offset], ldt, &work[work_offset], ldwork); 00548 00549 /* C := C - W * V */ 00550 00551 if (*n > *k) { 00552 00553 /* C2 := C2 - W * V2 */ 00554 00555 i__1 = *n - *k; 00556 template_blas_gemm("No transpose", "No transpose", m, &i__1, k, & 00557 c_b25, &work[work_offset], ldwork, &v_ref(1, *k + 00558 1), ldv, &c_b14, &c___ref(1, *k + 1), ldc); 00559 } 00560 00561 /* W := W * V1 */ 00562 00563 template_blas_trmm("Right", "Upper", "No transpose", "Unit", m, k, &c_b14, 00564 &v[v_offset], ldv, &work[work_offset], ldwork); 00565 00566 /* C1 := C1 - W */ 00567 00568 i__1 = *k; 00569 for (j = 1; j <= i__1; ++j) { 00570 i__2 = *m; 00571 for (i__ = 1; i__ <= i__2; ++i__) { 00572 c___ref(i__, j) = c___ref(i__, j) - work_ref(i__, j); 00573 /* L170: */ 00574 } 00575 /* L180: */ 00576 } 00577 00578 } 00579 00580 } else { 00581 00582 /* Let V = ( V1 V2 ) (V2: last K columns) 00583 where V2 is unit lower triangular. */ 00584 00585 if (template_blas_lsame(side, "L")) { 00586 00587 /* Form H * C or H' * C where C = ( C1 ) 00588 ( C2 ) 00589 00590 W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 00591 00592 W := C2' */ 00593 00594 i__1 = *k; 00595 for (j = 1; j <= i__1; ++j) { 00596 template_blas_copy(n, &c___ref(*m - *k + j, 1), ldc, &work_ref(1, j), 00597 &c__1); 00598 /* L190: */ 00599 } 00600 00601 /* W := W * V2' */ 00602 00603 template_blas_trmm("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, & 00604 v_ref(1, *m - *k + 1), ldv, &work[work_offset], 00605 ldwork); 00606 if (*m > *k) { 00607 00608 /* W := W + C1'*V1' */ 00609 00610 i__1 = *m - *k; 00611 template_blas_gemm("Transpose", "Transpose", n, k, &i__1, &c_b14, & 00612 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, & 00613 work[work_offset], ldwork); 00614 } 00615 00616 /* W := W * T' or W * T */ 00617 00618 template_blas_trmm("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[ 00619 t_offset], ldt, &work[work_offset], ldwork); 00620 00621 /* C := C - V' * W' */ 00622 00623 if (*m > *k) { 00624 00625 /* C1 := C1 - V1' * W' */ 00626 00627 i__1 = *m - *k; 00628 template_blas_gemm("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[ 00629 v_offset], ldv, &work[work_offset], ldwork, & 00630 c_b14, &c__[c_offset], ldc); 00631 } 00632 00633 /* W := W * V2 */ 00634 00635 template_blas_trmm("Right", "Lower", "No transpose", "Unit", n, k, &c_b14, 00636 &v_ref(1, *m - *k + 1), ldv, &work[work_offset], 00637 ldwork); 00638 00639 /* C2 := C2 - W' */ 00640 00641 i__1 = *k; 00642 for (j = 1; j <= i__1; ++j) { 00643 i__2 = *n; 00644 for (i__ = 1; i__ <= i__2; ++i__) { 00645 c___ref(*m - *k + j, i__) = c___ref(*m - *k + j, i__) 00646 - work_ref(i__, j); 00647 /* L200: */ 00648 } 00649 /* L210: */ 00650 } 00651 00652 } else if (template_blas_lsame(side, "R")) { 00653 00654 /* Form C * H or C * H' where C = ( C1 C2 ) 00655 00656 W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 00657 00658 W := C2 */ 00659 00660 i__1 = *k; 00661 for (j = 1; j <= i__1; ++j) { 00662 template_blas_copy(m, &c___ref(1, *n - *k + j), &c__1, &work_ref(1, j) 00663 , &c__1); 00664 /* L220: */ 00665 } 00666 00667 /* W := W * V2' */ 00668 00669 template_blas_trmm("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, & 00670 v_ref(1, *n - *k + 1), ldv, &work[work_offset], 00671 ldwork); 00672 if (*n > *k) { 00673 00674 /* W := W + C1 * V1' */ 00675 00676 i__1 = *n - *k; 00677 template_blas_gemm("No transpose", "Transpose", m, k, &i__1, &c_b14, & 00678 c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, & 00679 work[work_offset], ldwork); 00680 } 00681 00682 /* W := W * T or W * T' */ 00683 00684 template_blas_trmm("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[ 00685 t_offset], ldt, &work[work_offset], ldwork); 00686 00687 /* C := C - W * V */ 00688 00689 if (*n > *k) { 00690 00691 /* C1 := C1 - W * V1 */ 00692 00693 i__1 = *n - *k; 00694 template_blas_gemm("No transpose", "No transpose", m, &i__1, k, & 00695 c_b25, &work[work_offset], ldwork, &v[v_offset], 00696 ldv, &c_b14, &c__[c_offset], ldc); 00697 } 00698 00699 /* W := W * V2 */ 00700 00701 template_blas_trmm("Right", "Lower", "No transpose", "Unit", m, k, &c_b14, 00702 &v_ref(1, *n - *k + 1), ldv, &work[work_offset], 00703 ldwork); 00704 00705 /* C1 := C1 - W */ 00706 00707 i__1 = *k; 00708 for (j = 1; j <= i__1; ++j) { 00709 i__2 = *m; 00710 for (i__ = 1; i__ <= i__2; ++i__) { 00711 c___ref(i__, *n - *k + j) = c___ref(i__, *n - *k + j) 00712 - work_ref(i__, j); 00713 /* L230: */ 00714 } 00715 /* L240: */ 00716 } 00717 00718 } 00719 00720 } 00721 } 00722 00723 return 0; 00724 00725 /* End of DLARFB */ 00726 00727 } /* dlarfb_ */ 00728 00729 #undef v_ref 00730 #undef c___ref 00731 #undef work_ref 00732 00733 00734 #endif