00001 /* sspr.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int sspr_(char *uplo, integer *n, real *alpha, real *x, 00017 integer *incx, real *ap) 00018 { 00019 /* System generated locals */ 00020 integer i__1, i__2; 00021 00022 /* Local variables */ 00023 integer i__, j, k, kk, ix, jx, kx, info; 00024 real temp; 00025 extern logical lsame_(char *, char *); 00026 extern /* Subroutine */ int xerbla_(char *, integer *); 00027 00028 /* .. Scalar Arguments .. */ 00029 /* .. */ 00030 /* .. Array Arguments .. */ 00031 /* .. */ 00032 00033 /* Purpose */ 00034 /* ======= */ 00035 00036 /* SSPR performs the symmetric rank 1 operation */ 00037 00038 /* A := alpha*x*x' + A, */ 00039 00040 /* where alpha is a real scalar, x is an n element vector and A is an */ 00041 /* n by n symmetric matrix, supplied in packed form. */ 00042 00043 /* Arguments */ 00044 /* ========== */ 00045 00046 /* UPLO - CHARACTER*1. */ 00047 /* On entry, UPLO specifies whether the upper or lower */ 00048 /* triangular part of the matrix A is supplied in the packed */ 00049 /* array AP as follows: */ 00050 00051 /* UPLO = 'U' or 'u' The upper triangular part of A is */ 00052 /* supplied in AP. */ 00053 00054 /* UPLO = 'L' or 'l' The lower triangular part of A is */ 00055 /* supplied in AP. */ 00056 00057 /* Unchanged on exit. */ 00058 00059 /* N - INTEGER. */ 00060 /* On entry, N specifies the order of the matrix A. */ 00061 /* N must be at least zero. */ 00062 /* Unchanged on exit. */ 00063 00064 /* ALPHA - REAL . */ 00065 /* On entry, ALPHA specifies the scalar alpha. */ 00066 /* Unchanged on exit. */ 00067 00068 /* X - REAL array of dimension at least */ 00069 /* ( 1 + ( n - 1 )*abs( INCX ) ). */ 00070 /* Before entry, the incremented array X must contain the n */ 00071 /* element vector x. */ 00072 /* Unchanged on exit. */ 00073 00074 /* INCX - INTEGER. */ 00075 /* On entry, INCX specifies the increment for the elements of */ 00076 /* X. INCX must not be zero. */ 00077 /* Unchanged on exit. */ 00078 00079 /* AP - REAL array of DIMENSION at least */ 00080 /* ( ( n*( n + 1 ) )/2 ). */ 00081 /* Before entry with UPLO = 'U' or 'u', the array AP must */ 00082 /* contain the upper triangular part of the symmetric matrix */ 00083 /* packed sequentially, column by column, so that AP( 1 ) */ 00084 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */ 00085 /* and a( 2, 2 ) respectively, and so on. On exit, the array */ 00086 /* AP is overwritten by the upper triangular part of the */ 00087 /* updated matrix. */ 00088 /* Before entry with UPLO = 'L' or 'l', the array AP must */ 00089 /* contain the lower triangular part of the symmetric matrix */ 00090 /* packed sequentially, column by column, so that AP( 1 ) */ 00091 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */ 00092 /* and a( 3, 1 ) respectively, and so on. On exit, the array */ 00093 /* AP is overwritten by the lower triangular part of the */ 00094 /* updated matrix. */ 00095 00096 00097 /* Level 2 Blas routine. */ 00098 00099 /* -- Written on 22-October-1986. */ 00100 /* Jack Dongarra, Argonne National Lab. */ 00101 /* Jeremy Du Croz, Nag Central Office. */ 00102 /* Sven Hammarling, Nag Central Office. */ 00103 /* Richard Hanson, Sandia National Labs. */ 00104 00105 00106 /* .. Parameters .. */ 00107 /* .. */ 00108 /* .. Local Scalars .. */ 00109 /* .. */ 00110 /* .. External Functions .. */ 00111 /* .. */ 00112 /* .. External Subroutines .. */ 00113 /* .. */ 00114 00115 /* Test the input parameters. */ 00116 00117 /* Parameter adjustments */ 00118 --ap; 00119 --x; 00120 00121 /* Function Body */ 00122 info = 0; 00123 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) { 00124 info = 1; 00125 } else if (*n < 0) { 00126 info = 2; 00127 } else if (*incx == 0) { 00128 info = 5; 00129 } 00130 if (info != 0) { 00131 xerbla_("SSPR ", &info); 00132 return 0; 00133 } 00134 00135 /* Quick return if possible. */ 00136 00137 if (*n == 0 || *alpha == 0.f) { 00138 return 0; 00139 } 00140 00141 /* Set the start point in X if the increment is not unity. */ 00142 00143 if (*incx <= 0) { 00144 kx = 1 - (*n - 1) * *incx; 00145 } else if (*incx != 1) { 00146 kx = 1; 00147 } 00148 00149 /* Start the operations. In this version the elements of the array AP */ 00150 /* are accessed sequentially with one pass through AP. */ 00151 00152 kk = 1; 00153 if (lsame_(uplo, "U")) { 00154 00155 /* Form A when upper triangle is stored in AP. */ 00156 00157 if (*incx == 1) { 00158 i__1 = *n; 00159 for (j = 1; j <= i__1; ++j) { 00160 if (x[j] != 0.f) { 00161 temp = *alpha * x[j]; 00162 k = kk; 00163 i__2 = j; 00164 for (i__ = 1; i__ <= i__2; ++i__) { 00165 ap[k] += x[i__] * temp; 00166 ++k; 00167 /* L10: */ 00168 } 00169 } 00170 kk += j; 00171 /* L20: */ 00172 } 00173 } else { 00174 jx = kx; 00175 i__1 = *n; 00176 for (j = 1; j <= i__1; ++j) { 00177 if (x[jx] != 0.f) { 00178 temp = *alpha * x[jx]; 00179 ix = kx; 00180 i__2 = kk + j - 1; 00181 for (k = kk; k <= i__2; ++k) { 00182 ap[k] += x[ix] * temp; 00183 ix += *incx; 00184 /* L30: */ 00185 } 00186 } 00187 jx += *incx; 00188 kk += j; 00189 /* L40: */ 00190 } 00191 } 00192 } else { 00193 00194 /* Form A when lower triangle is stored in AP. */ 00195 00196 if (*incx == 1) { 00197 i__1 = *n; 00198 for (j = 1; j <= i__1; ++j) { 00199 if (x[j] != 0.f) { 00200 temp = *alpha * x[j]; 00201 k = kk; 00202 i__2 = *n; 00203 for (i__ = j; i__ <= i__2; ++i__) { 00204 ap[k] += x[i__] * temp; 00205 ++k; 00206 /* L50: */ 00207 } 00208 } 00209 kk = kk + *n - j + 1; 00210 /* L60: */ 00211 } 00212 } else { 00213 jx = kx; 00214 i__1 = *n; 00215 for (j = 1; j <= i__1; ++j) { 00216 if (x[jx] != 0.f) { 00217 temp = *alpha * x[jx]; 00218 ix = jx; 00219 i__2 = kk + *n - j; 00220 for (k = kk; k <= i__2; ++k) { 00221 ap[k] += x[ix] * temp; 00222 ix += *incx; 00223 /* L70: */ 00224 } 00225 } 00226 jx += *incx; 00227 kk = kk + *n - j + 1; 00228 /* L80: */ 00229 } 00230 } 00231 } 00232 00233 return 0; 00234 00235 /* End of SSPR . */ 00236 00237 } /* sspr_ */