template_blas_tpsv.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_BLAS_TPSV_HEADER
00036 #define TEMPLATE_BLAS_TPSV_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n, 
00041         const Treal *ap, Treal *x, const integer *incx)
00042 {
00043     /* System generated locals */
00044     integer i__1, i__2;
00045     /* Local variables */
00046      integer info;
00047      Treal temp;
00048      integer i__, j, k;
00049      integer kk, ix, jx, kx;
00050      logical nounit;
00051 /*  Purpose   
00052     =======   
00053     DTPSV  solves one of the systems of equations   
00054        A*x = b,   or   A'*x = b,   
00055     where b and x are n element vectors and A is an n by n unit, or   
00056     non-unit, upper or lower triangular matrix, supplied in packed form.   
00057     No test for singularity or near-singularity is included in this   
00058     routine. Such tests must be performed before calling this routine.   
00059     Parameters   
00060     ==========   
00061     UPLO   - CHARACTER*1.   
00062              On entry, UPLO specifies whether the matrix is an upper or   
00063              lower triangular matrix as follows:   
00064                 UPLO = 'U' or 'u'   A is an upper triangular matrix.   
00065                 UPLO = 'L' or 'l'   A is a lower triangular matrix.   
00066              Unchanged on exit.   
00067     TRANS  - CHARACTER*1.   
00068              On entry, TRANS specifies the equations to be solved as   
00069              follows:   
00070                 TRANS = 'N' or 'n'   A*x = b.   
00071                 TRANS = 'T' or 't'   A'*x = b.   
00072                 TRANS = 'C' or 'c'   A'*x = b.   
00073              Unchanged on exit.   
00074     DIAG   - CHARACTER*1.   
00075              On entry, DIAG specifies whether or not A is unit   
00076              triangular as follows:   
00077                 DIAG = 'U' or 'u'   A is assumed to be unit triangular.   
00078                 DIAG = 'N' or 'n'   A is not assumed to be unit   
00079                                     triangular.   
00080              Unchanged on exit.   
00081     N      - INTEGER.   
00082              On entry, N specifies the order of the matrix A.   
00083              N must be at least zero.   
00084              Unchanged on exit.   
00085     AP     - DOUBLE PRECISION array of DIMENSION at least   
00086              ( ( n*( n + 1 ) )/2 ).   
00087              Before entry with  UPLO = 'U' or 'u', the array AP must   
00088              contain the upper triangular matrix packed sequentially,   
00089              column by column, so that AP( 1 ) contains a( 1, 1 ),   
00090              AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )   
00091              respectively, and so on.   
00092              Before entry with UPLO = 'L' or 'l', the array AP must   
00093              contain the lower triangular matrix packed sequentially,   
00094              column by column, so that AP( 1 ) contains a( 1, 1 ),   
00095              AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )   
00096              respectively, and so on.   
00097              Note that when  DIAG = 'U' or 'u', the diagonal elements of   
00098              A are not referenced, but are assumed to be unity.   
00099              Unchanged on exit.   
00100     X      - DOUBLE PRECISION array of dimension at least   
00101              ( 1 + ( n - 1 )*abs( INCX ) ).   
00102              Before entry, the incremented array X must contain the n   
00103              element right-hand side vector b. On exit, X is overwritten   
00104              with the solution vector x.   
00105     INCX   - INTEGER.   
00106              On entry, INCX specifies the increment for the elements of   
00107              X. INCX must not be zero.   
00108              Unchanged on exit.   
00109     Level 2 Blas routine.   
00110     -- Written on 22-October-1986.   
00111        Jack Dongarra, Argonne National Lab.   
00112        Jeremy Du Croz, Nag Central Office.   
00113        Sven Hammarling, Nag Central Office.   
00114        Richard Hanson, Sandia National Labs.   
00115        Test the input parameters.   
00116        Parameter adjustments */
00117     --x;
00118     --ap;
00119     /* Initialization added by Elias to get rid of compiler warnings. */
00120     kx = 0;
00121     /* Function Body */
00122     info = 0;
00123     if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00124         info = 1;
00125     } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans, 
00126             "T") && ! template_blas_lsame(trans, "C")) {
00127         info = 2;
00128     } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag, 
00129             "N")) {
00130         info = 3;
00131     } else if (*n < 0) {
00132         info = 4;
00133     } else if (*incx == 0) {
00134         info = 7;
00135     }
00136     if (info != 0) {
00137         template_blas_erbla("DTPSV ", &info);
00138         return 0;
00139     }
00140 /*     Quick return if possible. */
00141     if (*n == 0) {
00142         return 0;
00143     }
00144     nounit = template_blas_lsame(diag, "N");
00145 /*     Set up the start point in X if the increment is not unity. This   
00146        will be  ( N - 1 )*INCX  too small for descending loops. */
00147     if (*incx <= 0) {
00148         kx = 1 - (*n - 1) * *incx;
00149     } else if (*incx != 1) {
00150         kx = 1;
00151     }
00152 /*     Start the operations. In this version the elements of AP are   
00153        accessed sequentially with one pass through AP. */
00154     if (template_blas_lsame(trans, "N")) {
00155 /*        Form  x := inv( A )*x. */
00156         if (template_blas_lsame(uplo, "U")) {
00157             kk = *n * (*n + 1) / 2;
00158             if (*incx == 1) {
00159                 for (j = *n; j >= 1; --j) {
00160                     if (x[j] != 0.) {
00161                         if (nounit) {
00162                             x[j] /= ap[kk];
00163                         }
00164                         temp = x[j];
00165                         k = kk - 1;
00166                         for (i__ = j - 1; i__ >= 1; --i__) {
00167                             x[i__] -= temp * ap[k];
00168                             --k;
00169 /* L10: */
00170                         }
00171                     }
00172                     kk -= j;
00173 /* L20: */
00174                 }
00175             } else {
00176                 jx = kx + (*n - 1) * *incx;
00177                 for (j = *n; j >= 1; --j) {
00178                     if (x[jx] != 0.) {
00179                         if (nounit) {
00180                             x[jx] /= ap[kk];
00181                         }
00182                         temp = x[jx];
00183                         ix = jx;
00184                         i__1 = kk - j + 1;
00185                         for (k = kk - 1; k >= i__1; --k) {
00186                             ix -= *incx;
00187                             x[ix] -= temp * ap[k];
00188 /* L30: */
00189                         }
00190                     }
00191                     jx -= *incx;
00192                     kk -= j;
00193 /* L40: */
00194                 }
00195             }
00196         } else {
00197             kk = 1;
00198             if (*incx == 1) {
00199                 i__1 = *n;
00200                 for (j = 1; j <= i__1; ++j) {
00201                     if (x[j] != 0.) {
00202                         if (nounit) {
00203                             x[j] /= ap[kk];
00204                         }
00205                         temp = x[j];
00206                         k = kk + 1;
00207                         i__2 = *n;
00208                         for (i__ = j + 1; i__ <= i__2; ++i__) {
00209                             x[i__] -= temp * ap[k];
00210                             ++k;
00211 /* L50: */
00212                         }
00213                     }
00214                     kk += *n - j + 1;
00215 /* L60: */
00216                 }
00217             } else {
00218                 jx = kx;
00219                 i__1 = *n;
00220                 for (j = 1; j <= i__1; ++j) {
00221                     if (x[jx] != 0.) {
00222                         if (nounit) {
00223                             x[jx] /= ap[kk];
00224                         }
00225                         temp = x[jx];
00226                         ix = jx;
00227                         i__2 = kk + *n - j;
00228                         for (k = kk + 1; k <= i__2; ++k) {
00229                             ix += *incx;
00230                             x[ix] -= temp * ap[k];
00231 /* L70: */
00232                         }
00233                     }
00234                     jx += *incx;
00235                     kk += *n - j + 1;
00236 /* L80: */
00237                 }
00238             }
00239         }
00240     } else {
00241 /*        Form  x := inv( A' )*x. */
00242         if (template_blas_lsame(uplo, "U")) {
00243             kk = 1;
00244             if (*incx == 1) {
00245                 i__1 = *n;
00246                 for (j = 1; j <= i__1; ++j) {
00247                     temp = x[j];
00248                     k = kk;
00249                     i__2 = j - 1;
00250                     for (i__ = 1; i__ <= i__2; ++i__) {
00251                         temp -= ap[k] * x[i__];
00252                         ++k;
00253 /* L90: */
00254                     }
00255                     if (nounit) {
00256                         temp /= ap[kk + j - 1];
00257                     }
00258                     x[j] = temp;
00259                     kk += j;
00260 /* L100: */
00261                 }
00262             } else {
00263                 jx = kx;
00264                 i__1 = *n;
00265                 for (j = 1; j <= i__1; ++j) {
00266                     temp = x[jx];
00267                     ix = kx;
00268                     i__2 = kk + j - 2;
00269                     for (k = kk; k <= i__2; ++k) {
00270                         temp -= ap[k] * x[ix];
00271                         ix += *incx;
00272 /* L110: */
00273                     }
00274                     if (nounit) {
00275                         temp /= ap[kk + j - 1];
00276                     }
00277                     x[jx] = temp;
00278                     jx += *incx;
00279                     kk += j;
00280 /* L120: */
00281                 }
00282             }
00283         } else {
00284             kk = *n * (*n + 1) / 2;
00285             if (*incx == 1) {
00286                 for (j = *n; j >= 1; --j) {
00287                     temp = x[j];
00288                     k = kk;
00289                     i__1 = j + 1;
00290                     for (i__ = *n; i__ >= i__1; --i__) {
00291                         temp -= ap[k] * x[i__];
00292                         --k;
00293 /* L130: */
00294                     }
00295                     if (nounit) {
00296                         temp /= ap[kk - *n + j];
00297                     }
00298                     x[j] = temp;
00299                     kk -= *n - j + 1;
00300 /* L140: */
00301                 }
00302             } else {
00303                 kx += (*n - 1) * *incx;
00304                 jx = kx;
00305                 for (j = *n; j >= 1; --j) {
00306                     temp = x[jx];
00307                     ix = kx;
00308                     i__1 = kk - (*n - (j + 1));
00309                     for (k = kk; k >= i__1; --k) {
00310                         temp -= ap[k] * x[ix];
00311                         ix -= *incx;
00312 /* L150: */
00313                     }
00314                     if (nounit) {
00315                         temp /= ap[kk - *n + j];
00316                     }
00317                     x[jx] = temp;
00318                     jx -= *incx;
00319                     kk -= *n - j + 1;
00320 /* L160: */
00321                 }
00322             }
00323         }
00324     }
00325     return 0;
00326 /*     End of DTPSV . */
00327 } /* dtpsv_ */
00328 
00329 #endif

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7