Logo ROOT   6.14/05
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 
15 namespace ROOT {
16 
17  namespace Minuit2 {
18 
19 
20 bool mnlsame(const char*, const char*);
21 int mnxerbla(const char*, int);
22 
23 int 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
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
Double_t x[n]
Definition: legend1.C:17
int mnxerbla(const char *, int)
Definition: mnxerbla.cxx:27
int mndspr(const char *, unsigned int, double, const double *, int, double *)
Definition: mndspr.cxx:23
bool mnlsame(const char *, const char *)
Definition: mnlsame.cxx:22
const Int_t n
Definition: legend1.C:16