00001 /* dspr2.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 dspr2_(char *uplo, integer *n, doublereal *alpha, 00017 doublereal *x, integer *incx, doublereal *y, integer *incy, 00018 doublereal *ap) 00019 { 00020 /* System generated locals */ 00021 integer i__1, i__2; 00022 00023 /* Local variables */ 00024 integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info; 00025 doublereal temp1, temp2; 00026 extern logical lsame_(char *, char *); 00027 extern /* Subroutine */ int xerbla_(char *, integer *); 00028 00029 /* .. Scalar Arguments .. */ 00030 /* .. */ 00031 /* .. Array Arguments .. */ 00032 /* .. */ 00033 00034 /* Purpose */ 00035 /* ======= */ 00036 00037 /* DSPR2 performs the symmetric rank 2 operation */ 00038 00039 /* A := alpha*x*y' + alpha*y*x' + A, */ 00040 00041 /* where alpha is a scalar, x and y are n element vectors and A is an */ 00042 /* n by n symmetric matrix, supplied in packed form. */ 00043 00044 /* Arguments */ 00045 /* ========== */ 00046 00047 /* UPLO - CHARACTER*1. */ 00048 /* On entry, UPLO specifies whether the upper or lower */ 00049 /* triangular part of the matrix A is supplied in the packed */ 00050 /* array AP as follows: */ 00051 00052 /* UPLO = 'U' or 'u' The upper triangular part of A is */ 00053 /* supplied in AP. */ 00054 00055 /* UPLO = 'L' or 'l' The lower triangular part of A is */ 00056 /* supplied in AP. */ 00057 00058 /* Unchanged on exit. */ 00059 00060 /* N - INTEGER. */ 00061 /* On entry, N specifies the order of the matrix A. */ 00062 /* N must be at least zero. */ 00063 /* Unchanged on exit. */ 00064 00065 /* ALPHA - DOUBLE PRECISION. */ 00066 /* On entry, ALPHA specifies the scalar alpha. */ 00067 /* Unchanged on exit. */ 00068 00069 /* X - DOUBLE PRECISION array of dimension at least */ 00070 /* ( 1 + ( n - 1 )*abs( INCX ) ). */ 00071 /* Before entry, the incremented array X must contain the n */ 00072 /* element vector x. */ 00073 /* Unchanged on exit. */ 00074 00075 /* INCX - INTEGER. */ 00076 /* On entry, INCX specifies the increment for the elements of */ 00077 /* X. INCX must not be zero. */ 00078 /* Unchanged on exit. */ 00079 00080 /* Y - DOUBLE PRECISION array of dimension at least */ 00081 /* ( 1 + ( n - 1 )*abs( INCY ) ). */ 00082 /* Before entry, the incremented array Y must contain the n */ 00083 /* element vector y. */ 00084 /* Unchanged on exit. */ 00085 00086 /* INCY - INTEGER. */ 00087 /* On entry, INCY specifies the increment for the elements of */ 00088 /* Y. INCY must not be zero. */ 00089 /* Unchanged on exit. */ 00090 00091 /* AP - DOUBLE PRECISION array of DIMENSION at least */ 00092 /* ( ( n*( n + 1 ) )/2 ). */ 00093 /* Before entry with UPLO = 'U' or 'u', the array AP must */ 00094 /* contain the upper triangular part of the symmetric matrix */ 00095 /* packed sequentially, column by column, so that AP( 1 ) */ 00096 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */ 00097 /* and a( 2, 2 ) respectively, and so on. On exit, the array */ 00098 /* AP is overwritten by the upper triangular part of the */ 00099 /* updated matrix. */ 00100 /* Before entry with UPLO = 'L' or 'l', the array AP must */ 00101 /* contain the lower triangular part of the symmetric matrix */ 00102 /* packed sequentially, column by column, so that AP( 1 ) */ 00103 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */ 00104 /* and a( 3, 1 ) respectively, and so on. On exit, the array */ 00105 /* AP is overwritten by the lower triangular part of the */ 00106 /* updated matrix. */ 00107 00108 00109 /* Level 2 Blas routine. */ 00110 00111 /* -- Written on 22-October-1986. */ 00112 /* Jack Dongarra, Argonne National Lab. */ 00113 /* Jeremy Du Croz, Nag Central Office. */ 00114 /* Sven Hammarling, Nag Central Office. */ 00115 /* Richard Hanson, Sandia National Labs. */ 00116 00117 00118 /* .. Parameters .. */ 00119 /* .. */ 00120 /* .. Local Scalars .. */ 00121 /* .. */ 00122 /* .. External Functions .. */ 00123 /* .. */ 00124 /* .. External Subroutines .. */ 00125 /* .. */ 00126 00127 /* Test the input parameters. */ 00128 00129 /* Parameter adjustments */ 00130 --ap; 00131 --y; 00132 --x; 00133 00134 /* Function Body */ 00135 info = 0; 00136 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) { 00137 info = 1; 00138 } else if (*n < 0) { 00139 info = 2; 00140 } else if (*incx == 0) { 00141 info = 5; 00142 } else if (*incy == 0) { 00143 info = 7; 00144 } 00145 if (info != 0) { 00146 xerbla_("DSPR2 ", &info); 00147 return 0; 00148 } 00149 00150 /* Quick return if possible. */ 00151 00152 if (*n == 0 || *alpha == 0.) { 00153 return 0; 00154 } 00155 00156 /* Set up the start points in X and Y if the increments are not both */ 00157 /* unity. */ 00158 00159 if (*incx != 1 || *incy != 1) { 00160 if (*incx > 0) { 00161 kx = 1; 00162 } else { 00163 kx = 1 - (*n - 1) * *incx; 00164 } 00165 if (*incy > 0) { 00166 ky = 1; 00167 } else { 00168 ky = 1 - (*n - 1) * *incy; 00169 } 00170 jx = kx; 00171 jy = ky; 00172 } 00173 00174 /* Start the operations. In this version the elements of the array AP */ 00175 /* are accessed sequentially with one pass through AP. */ 00176 00177 kk = 1; 00178 if (lsame_(uplo, "U")) { 00179 00180 /* Form A when upper triangle is stored in AP. */ 00181 00182 if (*incx == 1 && *incy == 1) { 00183 i__1 = *n; 00184 for (j = 1; j <= i__1; ++j) { 00185 if (x[j] != 0. || y[j] != 0.) { 00186 temp1 = *alpha * y[j]; 00187 temp2 = *alpha * x[j]; 00188 k = kk; 00189 i__2 = j; 00190 for (i__ = 1; i__ <= i__2; ++i__) { 00191 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2; 00192 ++k; 00193 /* L10: */ 00194 } 00195 } 00196 kk += j; 00197 /* L20: */ 00198 } 00199 } else { 00200 i__1 = *n; 00201 for (j = 1; j <= i__1; ++j) { 00202 if (x[jx] != 0. || y[jy] != 0.) { 00203 temp1 = *alpha * y[jy]; 00204 temp2 = *alpha * x[jx]; 00205 ix = kx; 00206 iy = ky; 00207 i__2 = kk + j - 1; 00208 for (k = kk; k <= i__2; ++k) { 00209 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2; 00210 ix += *incx; 00211 iy += *incy; 00212 /* L30: */ 00213 } 00214 } 00215 jx += *incx; 00216 jy += *incy; 00217 kk += j; 00218 /* L40: */ 00219 } 00220 } 00221 } else { 00222 00223 /* Form A when lower triangle is stored in AP. */ 00224 00225 if (*incx == 1 && *incy == 1) { 00226 i__1 = *n; 00227 for (j = 1; j <= i__1; ++j) { 00228 if (x[j] != 0. || y[j] != 0.) { 00229 temp1 = *alpha * y[j]; 00230 temp2 = *alpha * x[j]; 00231 k = kk; 00232 i__2 = *n; 00233 for (i__ = j; i__ <= i__2; ++i__) { 00234 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2; 00235 ++k; 00236 /* L50: */ 00237 } 00238 } 00239 kk = kk + *n - j + 1; 00240 /* L60: */ 00241 } 00242 } else { 00243 i__1 = *n; 00244 for (j = 1; j <= i__1; ++j) { 00245 if (x[jx] != 0. || y[jy] != 0.) { 00246 temp1 = *alpha * y[jy]; 00247 temp2 = *alpha * x[jx]; 00248 ix = jx; 00249 iy = jy; 00250 i__2 = kk + *n - j; 00251 for (k = kk; k <= i__2; ++k) { 00252 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2; 00253 ix += *incx; 00254 iy += *incy; 00255 /* L70: */ 00256 } 00257 } 00258 jx += *incx; 00259 jy += *incy; 00260 kk = kk + *n - j + 1; 00261 /* L80: */ 00262 } 00263 } 00264 } 00265 00266 return 0; 00267 00268 /* End of DSPR2 . */ 00269 00270 } /* dspr2_ */