Logo ROOT   6.16/01
Reference Guide
mndspr.cxx
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
10/* dspr.f -- translated by f2c (version 20010320).
11 You must link the resulting object file with the libraries:
12 -lf2c -lm (in that order)
13*/
14
15namespace ROOT {
16
17 namespace Minuit2 {
18
19
20bool mnlsame(const char*, const char*);
21int mnxerbla(const char*, int);
22
23int mndspr(const char* uplo, unsigned int n, double alpha,
24 const double* x, int incx, double* ap) {
25 /* System generated locals */
26 int i__1, i__2;
27
28 /* Local variables */
29 int info;
30 double temp;
31 int i__, j, k;
32 int kk, ix, jx, kx = 0;
33
34 /* .. Scalar Arguments .. */
35 /* .. Array Arguments .. */
36 /* .. */
37
38 /* Purpose */
39 /* ======= */
40
41 /* DSPR performs the symmetric rank 1 operation */
42
43 /* A := alpha*x*x' + A, */
44
45 /* where alpha is a real scalar, x is an n element vector and A is an */
46 /* n by n symmetric matrix, supplied in packed form. */
47
48 /* Parameters */
49 /* ========== */
50
51 /* UPLO - CHARACTER*1. */
52 /* On entry, UPLO specifies whether the Upper or Lower */
53 /* triangular part of the matrix A is supplied in the packed */
54 /* array AP as follows: */
55
56 /* UPLO = 'U' or 'u' The Upper triangular part of A is */
57 /* supplied in AP. */
58
59 /* UPLO = 'L' or 'l' The Lower triangular part of A is */
60 /* supplied in AP. */
61
62 /* Unchanged on exit. */
63
64 /* N - INTEGER. */
65 /* On entry, N specifies the order of the matrix A. */
66 /* N must be at least zero. */
67 /* Unchanged on exit. */
68
69 /* ALPHA - DOUBLE PRECISION. */
70 /* On entry, ALPHA specifies the scalar alpha. */
71 /* Unchanged on exit. */
72
73 /* X - DOUBLE PRECISION array of dimension at least */
74 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
75 /* Before entry, the incremented array X must contain the n */
76 /* element vector x. */
77 /* Unchanged on exit. */
78
79 /* INCX - INTEGER. */
80 /* On entry, INCX specifies the increment for the Elements of */
81 /* X. INCX must not be zero. */
82 /* Unchanged on exit. */
83
84 /* AP - DOUBLE PRECISION array of DIMENSION at least */
85 /* ( ( n*( n + 1 ) )/2 ). */
86 /* Before entry with UPLO = 'U' or 'u', the array AP must */
87 /* contain the Upper triangular part of the symmetric matrix */
88 /* packed sequentially, column by column, so that AP( 1 ) */
89 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
90 /* and a( 2, 2 ) respectively, and so on. On exit, the array */
91 /* AP is overwritten by the Upper triangular part of the */
92 /* updated matrix. */
93 /* Before entry with UPLO = 'L' or 'l', the array AP must */
94 /* contain the Lower triangular part of the symmetric matrix */
95 /* packed sequentially, column by column, so that AP( 1 ) */
96 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
97 /* and a( 3, 1 ) respectively, and so on. On exit, the array */
98 /* AP is overwritten by the Lower triangular part of the */
99 /* updated matrix. */
100
101
102 /* Level 2 Blas routine. */
103
104 /* -- Written on 22-October-1986. */
105 /* Jack Dongarra, Argonne National Lab. */
106 /* Jeremy Du Croz, Nag Central Office. */
107 /* Sven Hammarling, Nag Central Office. */
108 /* Richard Hanson, Sandia National Labs. */
109
110
111 /* .. Parameters .. */
112 /* .. Local Scalars .. */
113 /* .. External Functions .. */
114 /* .. External Subroutines .. */
115 /* .. */
116 /* .. Executable Statements .. */
117
118 /* Test the input parameters. */
119
120 /* Parameter adjustments */
121 --ap;
122 --x;
123
124 /* Function Body */
125 info = 0;
126 if (! mnlsame(uplo, "U") && ! mnlsame(uplo, "L")) {
127 info = 1;
128 }
129 // else if (n < 0) {
130 // info = 2;
131 // }
132 else if (incx == 0) {
133 info = 5;
134 }
135 if (info != 0) {
136 mnxerbla("DSPR ", info);
137 return 0;
138 }
139
140 /* Quick return if possible. */
141
142 if (n == 0 || alpha == 0.) {
143 return 0;
144 }
145
146 /* Set the start point in X if the increment is not unity. */
147
148 if (incx <= 0) {
149 kx = 1 - (n - 1) * incx;
150 } else if (incx != 1) {
151 kx = 1;
152 }
153
154 /* Start the operations. In this version the Elements of the array AP */
155 /* are accessed sequentially with one pass through AP. */
156
157 kk = 1;
158 if (mnlsame(uplo, "U")) {
159
160 /* Form A when Upper triangle is stored in AP. */
161
162 if (incx == 1) {
163 i__1 = n;
164 for (j = 1; j <= i__1; ++j) {
165 if (x[j] != 0.) {
166 temp = alpha * x[j];
167 k = kk;
168 i__2 = j;
169 for (i__ = 1; i__ <= i__2; ++i__) {
170 ap[k] += x[i__] * temp;
171 ++k;
172 /* L10: */
173 }
174 }
175 kk += j;
176 /* L20: */
177 }
178 } else {
179 jx = kx;
180 i__1 = n;
181 for (j = 1; j <= i__1; ++j) {
182 if (x[jx] != 0.) {
183 temp = alpha * x[jx];
184 ix = kx;
185 i__2 = kk + j - 1;
186 for (k = kk; k <= i__2; ++k) {
187 ap[k] += x[ix] * temp;
188 ix += incx;
189 /* L30: */
190 }
191 }
192 jx += incx;
193 kk += j;
194 /* L40: */
195 }
196 }
197 } else {
198
199 /* Form A when Lower triangle is stored in AP. */
200
201 if (incx == 1) {
202 i__1 = n;
203 for (j = 1; j <= i__1; ++j) {
204 if (x[j] != 0.) {
205 temp = alpha * x[j];
206 k = kk;
207 i__2 = n;
208 for (i__ = j; i__ <= i__2; ++i__) {
209 ap[k] += x[i__] * temp;
210 ++k;
211 /* L50: */
212 }
213 }
214 kk = kk + n - j + 1;
215 /* L60: */
216 }
217 } else {
218 jx = kx;
219 i__1 = n;
220 for (j = 1; j <= i__1; ++j) {
221 if (x[jx] != 0.) {
222 temp = alpha * x[jx];
223 ix = jx;
224 i__2 = kk + n - j;
225 for (k = kk; k <= i__2; ++k) {
226 ap[k] += x[ix] * temp;
227 ix += incx;
228 /* L70: */
229 }
230 }
231 jx += incx;
232 kk = kk + n - j + 1;
233 /* L80: */
234 }
235 }
236 }
237
238 return 0;
239
240 /* End of DSPR . */
241
242} /* dspr_ */
243
244
245 } // namespace Minuit2
246
247} // namespace ROOT
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
bool mnlsame(const char *, const char *)
Definition: mnlsame.cxx:22
int mnxerbla(const char *, int)
Definition: mnxerbla.cxx:27
int mndspr(const char *, unsigned int, double, const double *, int, double *)
Definition: mndspr.cxx:23
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21