| 1 | SUBROUTINE SGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, |
|---|
| 2 | $ WORK, LWORK, INFO ) |
|---|
| 3 | * |
|---|
| 4 | * -- LAPACK driver routine (version 3.3.1) -- |
|---|
| 5 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
|---|
| 6 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
|---|
| 7 | * -- April 2011 -- |
|---|
| 8 | * |
|---|
| 9 | * .. Scalar Arguments .. |
|---|
| 10 | INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
|---|
| 11 | REAL RCOND |
|---|
| 12 | * .. |
|---|
| 13 | * .. Array Arguments .. |
|---|
| 14 | INTEGER JPVT( * ) |
|---|
| 15 | REAL A( LDA, * ), B( LDB, * ), WORK( * ) |
|---|
| 16 | * .. |
|---|
| 17 | * |
|---|
| 18 | * Purpose |
|---|
| 19 | * ======= |
|---|
| 20 | * |
|---|
| 21 | * SGELSY computes the minimum-norm solution to a real linear least |
|---|
| 22 | * squares problem: |
|---|
| 23 | * minimize || A * X - B || |
|---|
| 24 | * using a complete orthogonal factorization of A. A is an M-by-N |
|---|
| 25 | * matrix which may be rank-deficient. |
|---|
| 26 | * |
|---|
| 27 | * Several right hand side vectors b and solution vectors x can be |
|---|
| 28 | * handled in a single call; they are stored as the columns of the |
|---|
| 29 | * M-by-NRHS right hand side matrix B and the N-by-NRHS solution |
|---|
| 30 | * matrix X. |
|---|
| 31 | * |
|---|
| 32 | * The routine first computes a QR factorization with column pivoting: |
|---|
| 33 | * A * P = Q * [ R11 R12 ] |
|---|
| 34 | * [ 0 R22 ] |
|---|
| 35 | * with R11 defined as the largest leading submatrix whose estimated |
|---|
| 36 | * condition number is less than 1/RCOND. The order of R11, RANK, |
|---|
| 37 | * is the effective rank of A. |
|---|
| 38 | * |
|---|
| 39 | * Then, R22 is considered to be negligible, and R12 is annihilated |
|---|
| 40 | * by orthogonal transformations from the right, arriving at the |
|---|
| 41 | * complete orthogonal factorization: |
|---|
| 42 | * A * P = Q * [ T11 0 ] * Z |
|---|
| 43 | * [ 0 0 ] |
|---|
| 44 | * The minimum-norm solution is then |
|---|
| 45 | * X = P * Z**T [ inv(T11)*Q1**T*B ] |
|---|
| 46 | * [ 0 ] |
|---|
| 47 | * where Q1 consists of the first RANK columns of Q. |
|---|
| 48 | * |
|---|
| 49 | * This routine is basically identical to the original xGELSX except |
|---|
| 50 | * three differences: |
|---|
| 51 | * o The call to the subroutine xGEQPF has been substituted by the |
|---|
| 52 | * the call to the subroutine xGEQP3. This subroutine is a Blas-3 |
|---|
| 53 | * version of the QR factorization with column pivoting. |
|---|
| 54 | * o Matrix B (the right hand side) is updated with Blas-3. |
|---|
| 55 | * o The permutation of matrix B (the right hand side) is faster and |
|---|
| 56 | * more simple. |
|---|
| 57 | * |
|---|
| 58 | * Arguments |
|---|
| 59 | * ========= |
|---|
| 60 | * |
|---|
| 61 | * M (input) INTEGER |
|---|
| 62 | * The number of rows of the matrix A. M >= 0. |
|---|
| 63 | * |
|---|
| 64 | * N (input) INTEGER |
|---|
| 65 | * The number of columns of the matrix A. N >= 0. |
|---|
| 66 | * |
|---|
| 67 | * NRHS (input) INTEGER |
|---|
| 68 | * The number of right hand sides, i.e., the number of |
|---|
| 69 | * columns of matrices B and X. NRHS >= 0. |
|---|
| 70 | * |
|---|
| 71 | * A (input/output) REAL array, dimension (LDA,N) |
|---|
| 72 | * On entry, the M-by-N matrix A. |
|---|
| 73 | * On exit, A has been overwritten by details of its |
|---|
| 74 | * complete orthogonal factorization. |
|---|
| 75 | * |
|---|
| 76 | * LDA (input) INTEGER |
|---|
| 77 | * The leading dimension of the array A. LDA >= max(1,M). |
|---|
| 78 | * |
|---|
| 79 | * B (input/output) REAL array, dimension (LDB,NRHS) |
|---|
| 80 | * On entry, the M-by-NRHS right hand side matrix B. |
|---|
| 81 | * On exit, the N-by-NRHS solution matrix X. |
|---|
| 82 | * |
|---|
| 83 | * LDB (input) INTEGER |
|---|
| 84 | * The leading dimension of the array B. LDB >= max(1,M,N). |
|---|
| 85 | * |
|---|
| 86 | * JPVT (input/output) INTEGER array, dimension (N) |
|---|
| 87 | * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted |
|---|
| 88 | * to the front of AP, otherwise column i is a free column. |
|---|
| 89 | * On exit, if JPVT(i) = k, then the i-th column of AP |
|---|
| 90 | * was the k-th column of A. |
|---|
| 91 | * |
|---|
| 92 | * RCOND (input) REAL |
|---|
| 93 | * RCOND is used to determine the effective rank of A, which |
|---|
| 94 | * is defined as the order of the largest leading triangular |
|---|
| 95 | * submatrix R11 in the QR factorization with pivoting of A, |
|---|
| 96 | * whose estimated condition number < 1/RCOND. |
|---|
| 97 | * |
|---|
| 98 | * RANK (output) INTEGER |
|---|
| 99 | * The effective rank of A, i.e., the order of the submatrix |
|---|
| 100 | * R11. This is the same as the order of the submatrix T11 |
|---|
| 101 | * in the complete orthogonal factorization of A. |
|---|
| 102 | * |
|---|
| 103 | * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) |
|---|
| 104 | * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|---|
| 105 | * |
|---|
| 106 | * LWORK (input) INTEGER |
|---|
| 107 | * The dimension of the array WORK. |
|---|
| 108 | * The unblocked strategy requires that: |
|---|
| 109 | * LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), |
|---|
| 110 | * where MN = min( M, N ). |
|---|
| 111 | * The block algorithm requires that: |
|---|
| 112 | * LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), |
|---|
| 113 | * where NB is an upper bound on the blocksize returned |
|---|
| 114 | * by ILAENV for the routines SGEQP3, STZRZF, STZRQF, SORMQR, |
|---|
| 115 | * and SORMRZ. |
|---|
| 116 | * |
|---|
| 117 | * If LWORK = -1, then a workspace query is assumed; the routine |
|---|
| 118 | * only calculates the optimal size of the WORK array, returns |
|---|
| 119 | * this value as the first entry of the WORK array, and no error |
|---|
| 120 | * message related to LWORK is issued by XERBLA. |
|---|
| 121 | * |
|---|
| 122 | * INFO (output) INTEGER |
|---|
| 123 | * = 0: successful exit |
|---|
| 124 | * < 0: If INFO = -i, the i-th argument had an illegal value. |
|---|
| 125 | * |
|---|
| 126 | * Further Details |
|---|
| 127 | * =============== |
|---|
| 128 | * |
|---|
| 129 | * Based on contributions by |
|---|
| 130 | * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA |
|---|
| 131 | * E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain |
|---|
| 132 | * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain |
|---|
| 133 | * |
|---|
| 134 | * ===================================================================== |
|---|
| 135 | * |
|---|
| 136 | * .. Parameters .. |
|---|
| 137 | INTEGER IMAX, IMIN |
|---|
| 138 | PARAMETER ( IMAX = 1, IMIN = 2 ) |
|---|
| 139 | REAL ZERO, ONE |
|---|
| 140 | PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) |
|---|
| 141 | * .. |
|---|
| 142 | * .. Local Scalars .. |
|---|
| 143 | LOGICAL LQUERY |
|---|
| 144 | INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN, |
|---|
| 145 | $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4 |
|---|
| 146 | REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, |
|---|
| 147 | $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE |
|---|
| 148 | * .. |
|---|
| 149 | * .. External Functions .. |
|---|
| 150 | INTEGER ILAENV |
|---|
| 151 | REAL SLAMCH, SLANGE |
|---|
| 152 | EXTERNAL ILAENV, SLAMCH, SLANGE |
|---|
| 153 | * .. |
|---|
| 154 | * .. External Subroutines .. |
|---|
| 155 | EXTERNAL SCOPY, SGEQP3, SLABAD, SLAIC1, SLASCL, SLASET, |
|---|
| 156 | $ SORMQR, SORMRZ, STRSM, STZRZF, XERBLA |
|---|
| 157 | * .. |
|---|
| 158 | * .. Intrinsic Functions .. |
|---|
| 159 | INTRINSIC ABS, MAX, MIN |
|---|
| 160 | * .. |
|---|
| 161 | * .. Executable Statements .. |
|---|
| 162 | * |
|---|
| 163 | MN = MIN( M, N ) |
|---|
| 164 | ISMIN = MN + 1 |
|---|
| 165 | ISMAX = 2*MN + 1 |
|---|
| 166 | * |
|---|
| 167 | * Test the input arguments. |
|---|
| 168 | * |
|---|
| 169 | INFO = 0 |
|---|
| 170 | LQUERY = ( LWORK.EQ.-1 ) |
|---|
| 171 | IF( M.LT.0 ) THEN |
|---|
| 172 | INFO = -1 |
|---|
| 173 | ELSE IF( N.LT.0 ) THEN |
|---|
| 174 | INFO = -2 |
|---|
| 175 | ELSE IF( NRHS.LT.0 ) THEN |
|---|
| 176 | INFO = -3 |
|---|
| 177 | ELSE IF( LDA.LT.MAX( 1, M ) ) THEN |
|---|
| 178 | INFO = -5 |
|---|
| 179 | ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN |
|---|
| 180 | INFO = -7 |
|---|
| 181 | END IF |
|---|
| 182 | * |
|---|
| 183 | * Figure out optimal block size |
|---|
| 184 | * |
|---|
| 185 | IF( INFO.EQ.0 ) THEN |
|---|
| 186 | IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN |
|---|
| 187 | LWKMIN = 1 |
|---|
| 188 | LWKOPT = 1 |
|---|
| 189 | ELSE |
|---|
| 190 | NB1 = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) |
|---|
| 191 | NB2 = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 ) |
|---|
| 192 | NB3 = ILAENV( 1, 'SORMQR', ' ', M, N, NRHS, -1 ) |
|---|
| 193 | NB4 = ILAENV( 1, 'SORMRQ', ' ', M, N, NRHS, -1 ) |
|---|
| 194 | NB = MAX( NB1, NB2, NB3, NB4 ) |
|---|
| 195 | LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS ) |
|---|
| 196 | LWKOPT = MAX( LWKMIN, |
|---|
| 197 | $ MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS ) |
|---|
| 198 | END IF |
|---|
| 199 | WORK( 1 ) = LWKOPT |
|---|
| 200 | * |
|---|
| 201 | IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN |
|---|
| 202 | INFO = -12 |
|---|
| 203 | END IF |
|---|
| 204 | END IF |
|---|
| 205 | * |
|---|
| 206 | IF( INFO.NE.0 ) THEN |
|---|
| 207 | CALL XERBLA( 'SGELSY', -INFO ) |
|---|
| 208 | RETURN |
|---|
| 209 | ELSE IF( LQUERY ) THEN |
|---|
| 210 | RETURN |
|---|
| 211 | END IF |
|---|
| 212 | * |
|---|
| 213 | * Quick return if possible |
|---|
| 214 | * |
|---|
| 215 | IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN |
|---|
| 216 | RANK = 0 |
|---|
| 217 | RETURN |
|---|
| 218 | END IF |
|---|
| 219 | * |
|---|
| 220 | * Get machine parameters |
|---|
| 221 | * |
|---|
| 222 | SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) |
|---|
| 223 | BIGNUM = ONE / SMLNUM |
|---|
| 224 | CALL SLABAD( SMLNUM, BIGNUM ) |
|---|
| 225 | * |
|---|
| 226 | * Scale A, B if max entries outside range [SMLNUM,BIGNUM] |
|---|
| 227 | * |
|---|
| 228 | ANRM = SLANGE( 'M', M, N, A, LDA, WORK ) |
|---|
| 229 | IASCL = 0 |
|---|
| 230 | IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN |
|---|
| 231 | * |
|---|
| 232 | * Scale matrix norm up to SMLNUM |
|---|
| 233 | * |
|---|
| 234 | CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) |
|---|
| 235 | IASCL = 1 |
|---|
| 236 | ELSE IF( ANRM.GT.BIGNUM ) THEN |
|---|
| 237 | * |
|---|
| 238 | * Scale matrix norm down to BIGNUM |
|---|
| 239 | * |
|---|
| 240 | CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) |
|---|
| 241 | IASCL = 2 |
|---|
| 242 | ELSE IF( ANRM.EQ.ZERO ) THEN |
|---|
| 243 | * |
|---|
| 244 | * Matrix all zero. Return zero solution. |
|---|
| 245 | * |
|---|
| 246 | CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) |
|---|
| 247 | RANK = 0 |
|---|
| 248 | GO TO 70 |
|---|
| 249 | END IF |
|---|
| 250 | * |
|---|
| 251 | BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK ) |
|---|
| 252 | IBSCL = 0 |
|---|
| 253 | IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN |
|---|
| 254 | * |
|---|
| 255 | * Scale matrix norm up to SMLNUM |
|---|
| 256 | * |
|---|
| 257 | CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) |
|---|
| 258 | IBSCL = 1 |
|---|
| 259 | ELSE IF( BNRM.GT.BIGNUM ) THEN |
|---|
| 260 | * |
|---|
| 261 | * Scale matrix norm down to BIGNUM |
|---|
| 262 | * |
|---|
| 263 | CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) |
|---|
| 264 | IBSCL = 2 |
|---|
| 265 | END IF |
|---|
| 266 | * |
|---|
| 267 | * Compute QR factorization with column pivoting of A: |
|---|
| 268 | * A * P = Q * R |
|---|
| 269 | * |
|---|
| 270 | CALL SGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), |
|---|
| 271 | $ LWORK-MN, INFO ) |
|---|
| 272 | WSIZE = MN + WORK( MN+1 ) |
|---|
| 273 | * |
|---|
| 274 | * workspace: MN+2*N+NB*(N+1). |
|---|
| 275 | * Details of Householder rotations stored in WORK(1:MN). |
|---|
| 276 | * |
|---|
| 277 | * Determine RANK using incremental condition estimation |
|---|
| 278 | * |
|---|
| 279 | WORK( ISMIN ) = ONE |
|---|
| 280 | WORK( ISMAX ) = ONE |
|---|
| 281 | SMAX = ABS( A( 1, 1 ) ) |
|---|
| 282 | SMIN = SMAX |
|---|
| 283 | IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN |
|---|
| 284 | RANK = 0 |
|---|
| 285 | CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) |
|---|
| 286 | GO TO 70 |
|---|
| 287 | ELSE |
|---|
| 288 | RANK = 1 |
|---|
| 289 | END IF |
|---|
| 290 | * |
|---|
| 291 | 10 CONTINUE |
|---|
| 292 | IF( RANK.LT.MN ) THEN |
|---|
| 293 | I = RANK + 1 |
|---|
| 294 | CALL SLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), |
|---|
| 295 | $ A( I, I ), SMINPR, S1, C1 ) |
|---|
| 296 | CALL SLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), |
|---|
| 297 | $ A( I, I ), SMAXPR, S2, C2 ) |
|---|
| 298 | * |
|---|
| 299 | IF( SMAXPR*RCOND.LE.SMINPR ) THEN |
|---|
| 300 | DO 20 I = 1, RANK |
|---|
| 301 | WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) |
|---|
| 302 | WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) |
|---|
| 303 | 20 CONTINUE |
|---|
| 304 | WORK( ISMIN+RANK ) = C1 |
|---|
| 305 | WORK( ISMAX+RANK ) = C2 |
|---|
| 306 | SMIN = SMINPR |
|---|
| 307 | SMAX = SMAXPR |
|---|
| 308 | RANK = RANK + 1 |
|---|
| 309 | GO TO 10 |
|---|
| 310 | END IF |
|---|
| 311 | END IF |
|---|
| 312 | * |
|---|
| 313 | * workspace: 3*MN. |
|---|
| 314 | * |
|---|
| 315 | * Logically partition R = [ R11 R12 ] |
|---|
| 316 | * [ 0 R22 ] |
|---|
| 317 | * where R11 = R(1:RANK,1:RANK) |
|---|
| 318 | * |
|---|
| 319 | * [R11,R12] = [ T11, 0 ] * Y |
|---|
| 320 | * |
|---|
| 321 | IF( RANK.LT.N ) |
|---|
| 322 | $ CALL STZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ), |
|---|
| 323 | $ LWORK-2*MN, INFO ) |
|---|
| 324 | * |
|---|
| 325 | * workspace: 2*MN. |
|---|
| 326 | * Details of Householder rotations stored in WORK(MN+1:2*MN) |
|---|
| 327 | * |
|---|
| 328 | * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) |
|---|
| 329 | * |
|---|
| 330 | CALL SORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), |
|---|
| 331 | $ B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) |
|---|
| 332 | WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) ) |
|---|
| 333 | * |
|---|
| 334 | * workspace: 2*MN+NB*NRHS. |
|---|
| 335 | * |
|---|
| 336 | * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) |
|---|
| 337 | * |
|---|
| 338 | CALL STRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, |
|---|
| 339 | $ NRHS, ONE, A, LDA, B, LDB ) |
|---|
| 340 | * |
|---|
| 341 | DO 40 J = 1, NRHS |
|---|
| 342 | DO 30 I = RANK + 1, N |
|---|
| 343 | B( I, J ) = ZERO |
|---|
| 344 | 30 CONTINUE |
|---|
| 345 | 40 CONTINUE |
|---|
| 346 | * |
|---|
| 347 | * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS) |
|---|
| 348 | * |
|---|
| 349 | IF( RANK.LT.N ) THEN |
|---|
| 350 | CALL SORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A, |
|---|
| 351 | $ LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ), |
|---|
| 352 | $ LWORK-2*MN, INFO ) |
|---|
| 353 | END IF |
|---|
| 354 | * |
|---|
| 355 | * workspace: 2*MN+NRHS. |
|---|
| 356 | * |
|---|
| 357 | * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) |
|---|
| 358 | * |
|---|
| 359 | DO 60 J = 1, NRHS |
|---|
| 360 | DO 50 I = 1, N |
|---|
| 361 | WORK( JPVT( I ) ) = B( I, J ) |
|---|
| 362 | 50 CONTINUE |
|---|
| 363 | CALL SCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 ) |
|---|
| 364 | 60 CONTINUE |
|---|
| 365 | * |
|---|
| 366 | * workspace: N. |
|---|
| 367 | * |
|---|
| 368 | * Undo scaling |
|---|
| 369 | * |
|---|
| 370 | IF( IASCL.EQ.1 ) THEN |
|---|
| 371 | CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) |
|---|
| 372 | CALL SLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, |
|---|
| 373 | $ INFO ) |
|---|
| 374 | ELSE IF( IASCL.EQ.2 ) THEN |
|---|
| 375 | CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) |
|---|
| 376 | CALL SLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, |
|---|
| 377 | $ INFO ) |
|---|
| 378 | END IF |
|---|
| 379 | IF( IBSCL.EQ.1 ) THEN |
|---|
| 380 | CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) |
|---|
| 381 | ELSE IF( IBSCL.EQ.2 ) THEN |
|---|
| 382 | CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) |
|---|
| 383 | END IF |
|---|
| 384 | * |
|---|
| 385 | 70 CONTINUE |
|---|
| 386 | WORK( 1 ) = LWKOPT |
|---|
| 387 | * |
|---|
| 388 | RETURN |
|---|
| 389 | * |
|---|
| 390 | * End of SGELSY |
|---|
| 391 | * |
|---|
| 392 | END |
|---|
![(please configure the [header_logo] section in trac.ini)](/grisli/chrome/site/Laurentide-16km-3D-Umag-water-bordure-GRISLI-small.png)
