template_lapack_lasr.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_LASR_HEADER
00036 #define TEMPLATE_LAPACK_LASR_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m,
00041          const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer *
00042         lda)
00043 {
00044 /*  -- LAPACK auxiliary routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        October 31, 1992   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DLASR   performs the transformation   
00054 
00055        A := P*A,   when SIDE = 'L' or 'l'  (  Left-hand side )   
00056 
00057        A := A*P',  when SIDE = 'R' or 'r'  ( Right-hand side )   
00058 
00059     where A is an m by n real matrix and P is an orthogonal matrix,   
00060     consisting of a sequence of plane rotations determined by the   
00061     parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'   
00062     and z = n when SIDE = 'R' or 'r' ):   
00063 
00064     When  DIRECT = 'F' or 'f'  ( Forward sequence ) then   
00065 
00066        P = P( z - 1 )*...*P( 2 )*P( 1 ),   
00067 
00068     and when DIRECT = 'B' or 'b'  ( Backward sequence ) then   
00069 
00070        P = P( 1 )*P( 2 )*...*P( z - 1 ),   
00071 
00072     where  P( k ) is a plane rotation matrix for the following planes:   
00073 
00074        when  PIVOT = 'V' or 'v'  ( Variable pivot ),   
00075           the plane ( k, k + 1 )   
00076 
00077        when  PIVOT = 'T' or 't'  ( Top pivot ),   
00078           the plane ( 1, k + 1 )   
00079 
00080        when  PIVOT = 'B' or 'b'  ( Bottom pivot ),   
00081           the plane ( k, z )   
00082 
00083     c( k ) and s( k )  must contain the  cosine and sine that define the   
00084     matrix  P( k ).  The two by two plane rotation part of the matrix   
00085     P( k ), R( k ), is assumed to be of the form   
00086 
00087        R( k ) = (  c( k )  s( k ) ).   
00088                 ( -s( k )  c( k ) )   
00089 
00090     This version vectorises across rows of the array A when SIDE = 'L'.   
00091 
00092     Arguments   
00093     =========   
00094 
00095     SIDE    (input) CHARACTER*1   
00096             Specifies whether the plane rotation matrix P is applied to   
00097             A on the left or the right.   
00098             = 'L':  Left, compute A := P*A   
00099             = 'R':  Right, compute A:= A*P'   
00100 
00101     DIRECT  (input) CHARACTER*1   
00102             Specifies whether P is a forward or backward sequence of   
00103             plane rotations.   
00104             = 'F':  Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )   
00105             = 'B':  Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )   
00106 
00107     PIVOT   (input) CHARACTER*1   
00108             Specifies the plane for which P(k) is a plane rotation   
00109             matrix.   
00110             = 'V':  Variable pivot, the plane (k,k+1)   
00111             = 'T':  Top pivot, the plane (1,k+1)   
00112             = 'B':  Bottom pivot, the plane (k,z)   
00113 
00114     M       (input) INTEGER   
00115             The number of rows of the matrix A.  If m <= 1, an immediate   
00116             return is effected.   
00117 
00118     N       (input) INTEGER   
00119             The number of columns of the matrix A.  If n <= 1, an   
00120             immediate return is effected.   
00121 
00122     C, S    (input) DOUBLE PRECISION arrays, dimension   
00123                     (M-1) if SIDE = 'L'   
00124                     (N-1) if SIDE = 'R'   
00125             c(k) and s(k) contain the cosine and sine that define the   
00126             matrix P(k).  The two by two plane rotation part of the   
00127             matrix P(k), R(k), is assumed to be of the form   
00128             R( k ) = (  c( k )  s( k ) ).   
00129                      ( -s( k )  c( k ) )   
00130 
00131     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00132             The m by n matrix A.  On exit, A is overwritten by P*A if   
00133             SIDE = 'R' or by A*P' if SIDE = 'L'.   
00134 
00135     LDA     (input) INTEGER   
00136             The leading dimension of the array A.  LDA >= max(1,M).   
00137 
00138     =====================================================================   
00139 
00140 
00141        Test the input parameters   
00142 
00143        Parameter adjustments */
00144     /* System generated locals */
00145     integer a_dim1, a_offset, i__1, i__2;
00146     /* Local variables */
00147      integer info;
00148      Treal temp;
00149      integer i__, j;
00150      Treal ctemp, stemp;
00151 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00152 
00153     --c__;
00154     --s;
00155     a_dim1 = *lda;
00156     a_offset = 1 + a_dim1 * 1;
00157     a -= a_offset;
00158 
00159     /* Function Body */
00160     info = 0;
00161     if (! (template_blas_lsame(side, "L") || template_blas_lsame(side, "R"))) {
00162         info = 1;
00163     } else if (! (template_blas_lsame(pivot, "V") || template_blas_lsame(pivot, 
00164             "T") || template_blas_lsame(pivot, "B"))) {
00165         info = 2;
00166     } else if (! (template_blas_lsame(direct, "F") || template_blas_lsame(direct, 
00167             "B"))) {
00168         info = 3;
00169     } else if (*m < 0) {
00170         info = 4;
00171     } else if (*n < 0) {
00172         info = 5;
00173     } else if (*lda < maxMACRO(1,*m)) {
00174         info = 9;
00175     }
00176     if (info != 0) {
00177         template_blas_erbla("LASR  ", &info);
00178         return 0;
00179     }
00180 
00181 /*     Quick return if possible */
00182 
00183     if (*m == 0 || *n == 0) {
00184         return 0;
00185     }
00186     if (template_blas_lsame(side, "L")) {
00187 
00188 /*        Form  P * A */
00189 
00190         if (template_blas_lsame(pivot, "V")) {
00191             if (template_blas_lsame(direct, "F")) {
00192                 i__1 = *m - 1;
00193                 for (j = 1; j <= i__1; ++j) {
00194                     ctemp = c__[j];
00195                     stemp = s[j];
00196                     if (ctemp != 1. || stemp != 0.) {
00197                         i__2 = *n;
00198                         for (i__ = 1; i__ <= i__2; ++i__) {
00199                             temp = a_ref(j + 1, i__);
00200                             a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
00201                                     j, i__);
00202                             a_ref(j, i__) = stemp * temp + ctemp * a_ref(j, 
00203                                     i__);
00204 /* L10: */
00205                         }
00206                     }
00207 /* L20: */
00208                 }
00209             } else if (template_blas_lsame(direct, "B")) {
00210                 for (j = *m - 1; j >= 1; --j) {
00211                     ctemp = c__[j];
00212                     stemp = s[j];
00213                     if (ctemp != 1. || stemp != 0.) {
00214                         i__1 = *n;
00215                         for (i__ = 1; i__ <= i__1; ++i__) {
00216                             temp = a_ref(j + 1, i__);
00217                             a_ref(j + 1, i__) = ctemp * temp - stemp * a_ref(
00218                                     j, i__);
00219                             a_ref(j, i__) = stemp * temp + ctemp * a_ref(j, 
00220                                     i__);
00221 /* L30: */
00222                         }
00223                     }
00224 /* L40: */
00225                 }
00226             }
00227         } else if (template_blas_lsame(pivot, "T")) {
00228             if (template_blas_lsame(direct, "F")) {
00229                 i__1 = *m;
00230                 for (j = 2; j <= i__1; ++j) {
00231                     ctemp = c__[j - 1];
00232                     stemp = s[j - 1];
00233                     if (ctemp != 1. || stemp != 0.) {
00234                         i__2 = *n;
00235                         for (i__ = 1; i__ <= i__2; ++i__) {
00236                             temp = a_ref(j, i__);
00237                             a_ref(j, i__) = ctemp * temp - stemp * a_ref(1, 
00238                                     i__);
00239                             a_ref(1, i__) = stemp * temp + ctemp * a_ref(1, 
00240                                     i__);
00241 /* L50: */
00242                         }
00243                     }
00244 /* L60: */
00245                 }
00246             } else if (template_blas_lsame(direct, "B")) {
00247                 for (j = *m; j >= 2; --j) {
00248                     ctemp = c__[j - 1];
00249                     stemp = s[j - 1];
00250                     if (ctemp != 1. || stemp != 0.) {
00251                         i__1 = *n;
00252                         for (i__ = 1; i__ <= i__1; ++i__) {
00253                             temp = a_ref(j, i__);
00254                             a_ref(j, i__) = ctemp * temp - stemp * a_ref(1, 
00255                                     i__);
00256                             a_ref(1, i__) = stemp * temp + ctemp * a_ref(1, 
00257                                     i__);
00258 /* L70: */
00259                         }
00260                     }
00261 /* L80: */
00262                 }
00263             }
00264         } else if (template_blas_lsame(pivot, "B")) {
00265             if (template_blas_lsame(direct, "F")) {
00266                 i__1 = *m - 1;
00267                 for (j = 1; j <= i__1; ++j) {
00268                     ctemp = c__[j];
00269                     stemp = s[j];
00270                     if (ctemp != 1. || stemp != 0.) {
00271                         i__2 = *n;
00272                         for (i__ = 1; i__ <= i__2; ++i__) {
00273                             temp = a_ref(j, i__);
00274                             a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp * 
00275                                     temp;
00276                             a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp * 
00277                                     temp;
00278 /* L90: */
00279                         }
00280                     }
00281 /* L100: */
00282                 }
00283             } else if (template_blas_lsame(direct, "B")) {
00284                 for (j = *m - 1; j >= 1; --j) {
00285                     ctemp = c__[j];
00286                     stemp = s[j];
00287                     if (ctemp != 1. || stemp != 0.) {
00288                         i__1 = *n;
00289                         for (i__ = 1; i__ <= i__1; ++i__) {
00290                             temp = a_ref(j, i__);
00291                             a_ref(j, i__) = stemp * a_ref(*m, i__) + ctemp * 
00292                                     temp;
00293                             a_ref(*m, i__) = ctemp * a_ref(*m, i__) - stemp * 
00294                                     temp;
00295 /* L110: */
00296                         }
00297                     }
00298 /* L120: */
00299                 }
00300             }
00301         }
00302     } else if (template_blas_lsame(side, "R")) {
00303 
00304 /*        Form A * P' */
00305 
00306         if (template_blas_lsame(pivot, "V")) {
00307             if (template_blas_lsame(direct, "F")) {
00308                 i__1 = *n - 1;
00309                 for (j = 1; j <= i__1; ++j) {
00310                     ctemp = c__[j];
00311                     stemp = s[j];
00312                     if (ctemp != 1. || stemp != 0.) {
00313                         i__2 = *m;
00314                         for (i__ = 1; i__ <= i__2; ++i__) {
00315                             temp = a_ref(i__, j + 1);
00316                             a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
00317                                     i__, j);
00318                             a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__, 
00319                                     j);
00320 /* L130: */
00321                         }
00322                     }
00323 /* L140: */
00324                 }
00325             } else if (template_blas_lsame(direct, "B")) {
00326                 for (j = *n - 1; j >= 1; --j) {
00327                     ctemp = c__[j];
00328                     stemp = s[j];
00329                     if (ctemp != 1. || stemp != 0.) {
00330                         i__1 = *m;
00331                         for (i__ = 1; i__ <= i__1; ++i__) {
00332                             temp = a_ref(i__, j + 1);
00333                             a_ref(i__, j + 1) = ctemp * temp - stemp * a_ref(
00334                                     i__, j);
00335                             a_ref(i__, j) = stemp * temp + ctemp * a_ref(i__, 
00336                                     j);
00337 /* L150: */
00338                         }
00339                     }
00340 /* L160: */
00341                 }
00342             }
00343         } else if (template_blas_lsame(pivot, "T")) {
00344             if (template_blas_lsame(direct, "F")) {
00345                 i__1 = *n;
00346                 for (j = 2; j <= i__1; ++j) {
00347                     ctemp = c__[j - 1];
00348                     stemp = s[j - 1];
00349                     if (ctemp != 1. || stemp != 0.) {
00350                         i__2 = *m;
00351                         for (i__ = 1; i__ <= i__2; ++i__) {
00352                             temp = a_ref(i__, j);
00353                             a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__, 
00354                                     1);
00355                             a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__, 
00356                                     1);
00357 /* L170: */
00358                         }
00359                     }
00360 /* L180: */
00361                 }
00362             } else if (template_blas_lsame(direct, "B")) {
00363                 for (j = *n; j >= 2; --j) {
00364                     ctemp = c__[j - 1];
00365                     stemp = s[j - 1];
00366                     if (ctemp != 1. || stemp != 0.) {
00367                         i__1 = *m;
00368                         for (i__ = 1; i__ <= i__1; ++i__) {
00369                             temp = a_ref(i__, j);
00370                             a_ref(i__, j) = ctemp * temp - stemp * a_ref(i__, 
00371                                     1);
00372                             a_ref(i__, 1) = stemp * temp + ctemp * a_ref(i__, 
00373                                     1);
00374 /* L190: */
00375                         }
00376                     }
00377 /* L200: */
00378                 }
00379             }
00380         } else if (template_blas_lsame(pivot, "B")) {
00381             if (template_blas_lsame(direct, "F")) {
00382                 i__1 = *n - 1;
00383                 for (j = 1; j <= i__1; ++j) {
00384                     ctemp = c__[j];
00385                     stemp = s[j];
00386                     if (ctemp != 1. || stemp != 0.) {
00387                         i__2 = *m;
00388                         for (i__ = 1; i__ <= i__2; ++i__) {
00389                             temp = a_ref(i__, j);
00390                             a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp * 
00391                                     temp;
00392                             a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp * 
00393                                     temp;
00394 /* L210: */
00395                         }
00396                     }
00397 /* L220: */
00398                 }
00399             } else if (template_blas_lsame(direct, "B")) {
00400                 for (j = *n - 1; j >= 1; --j) {
00401                     ctemp = c__[j];
00402                     stemp = s[j];
00403                     if (ctemp != 1. || stemp != 0.) {
00404                         i__1 = *m;
00405                         for (i__ = 1; i__ <= i__1; ++i__) {
00406                             temp = a_ref(i__, j);
00407                             a_ref(i__, j) = stemp * a_ref(i__, *n) + ctemp * 
00408                                     temp;
00409                             a_ref(i__, *n) = ctemp * a_ref(i__, *n) - stemp * 
00410                                     temp;
00411 /* L230: */
00412                         }
00413                     }
00414 /* L240: */
00415                 }
00416             }
00417         }
00418     }
00419 
00420     return 0;
00421 
00422 /*     End of DLASR */
00423 
00424 } /* dlasr_ */
00425 
00426 #undef a_ref
00427 
00428 
00429 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1