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