template_lapack_larfb.h

Go to the documentation of this file.
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

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7