ztpgqrt.c 15.6 KB
Newer Older
1
/**
2 3
 *
 * @file ztpgqrt.c
4
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
5 6
 * @copyright 2009-2016 The University of Tennessee and The University of
 *                      Tennessee Research Foundation. All rights reserved.
7 8 9
 * @copyright 2012-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 * @copyright 2016-2018 KAUST. All rights reserved.
10
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
11
 ***
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
 *
 *
 *  MORSE computational routines
 *  MORSE is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 0.9.0
 * @author Mathieu Faverge
 * @date 2016-12-15
 * @precisions normal z -> s d c
 *
 **/
#include "control/common.h"

/**
 ******************************************************************************
 *
 * @ingroup MORSE_Complex64_t
 *
 *  MORSE_ztpgqrt - Generates a partial Q matrix formed with a blocked QR
 *  factorization of a "triangular-pentagonal" matrix C, which is composed of an
 *  unused triangular block and a pentagonal block V, using the compact
 *  representation for Q. See MORSE_ztpqrt() to generate V.
 *
 *******************************************************************************
 *
 * @param[in] M
39
 *          The number of rows of the matrices Q2, and V2. M >= 0.
40 41
 *
 * @param[in] N
42
 *          The number of columns of the matrices Q1, and Q1. N >= 0.
43 44 45 46
 *
 * @param[in] K
 *          The number of elementary reflectors whose product defines
 *          the matrix Q in the matrix V.
47 48
 *          The number of rows of the matrices V1, and Q1.
 *          The number of columns of the matrices V1, and V2.
49 50
 *
 * @param[in] L
51
 *          The number of rows of the upper trapezoidal part of V2.
52 53
 *          MIN(M,N) >= L >= 0.  See Further Details.
 *
54
 * @param[in] V1
55 56
 *          The i-th row must contain the vector which defines the
 *          elementary reflector H(i), for i = 1,2,...,k, as returned by
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
 *          MORSE_ztpqrt().
 *          V1 is a matrix of size K-by-K.
 *
 * @param[in] LDV1
 *          The leading dimension of the array V1. LDV1 >= max(1,K).
 *
 * @param[int] descT1
 *          The auxiliary factorization data generated by the call to
 *          MORSE_zgeqrf() on V1.
 *
 * @param[in] V2
 *          The i-th row must contain the vector which defines the
 *          elementary reflector H(i), for i = 1,2,...,k, as returned by
 *          MORSE_ztpqrt() in the first k rows of its array argument V2.
 *          V2 is a matrix of size M-by-K. The first M-L rows
72 73
 *          are rectangular, and the last L rows are upper trapezoidal.
 *
74 75
 * @param[in] LDV2
 *          The leading dimension of the array V2. LDV2 >= max(1,M).
76
 *
77 78 79
 * @param[int] descT2
 *          The auxiliary factorization data, generated by the call to
 *          MORSE_ztpqrt() on V1 and V2, and associated to V2.
80
 *
81 82 83 84 85
 * @param[in,out] Q1
 *          Q1 is COMPLEX*16 array, dimension (LDQ1,N)
 *          On entry, the K-by-N matrix Q1.
 *          On exit, Q1 is overwritten by the corresponding block of
 *          Q*Q1.  See Further Details.
86
 *
87 88
 * @param[in] LDQ1
 *          The leading dimension of the array Q1. LDQ1 >= max(1,K).
89
 *
90 91 92
 * @param[in,out] Q2
 *          On entry, the pentagonal M-by-N matrix Q2.
 *          On exit, Q2 contains Q.
93
 *
94 95
 * @param[in] LDQ2
 *          The leading dimension of the array Q2.  LDQ2 >= max(1,M).
96 97 98 99 100 101
 *
 * @par Further Details:
 * =====================
 *
 *  The input matrix Q is a (K+M)-by-N matrix
 *
102 103
 *               Q = [ Q1 ]
 *                   [ Q2 ]
104
 *
105 106 107 108
 *  where Q1 is an identity matrix, and Q2 is a M-by-N matrix of 0.
 *  V is a matrix of householder reflectors with a pentagonal shape consisting
 *  of a K-by-K rectangular matrix V1 on top of a matrix V2 composed of
 *  (M-L)-by-K rectangular part on top of a L-by-N upper trapezoidal matrix:
109
 *
110 111 112
 *               V = [ V1  ]  <-     K-by-K rectangular
 *                   [ V2a ]  <- (M-L)-by-K rectangular
 *                   [ V2b ]  <-     L-by-K upper trapezoidal.
113
 *
114 115 116 117
 *  The upper trapezoidal part of the matrix V2 consists of the first L rows of
 *  a K-by-K upper triangular matrix, where 0 <= L <= MIN(M,K).  If L=0, V2 is
 *  rectangular M-by-K; if M=L=K, V2 is upper triangular. Those are the two
 *  cases only handled for now.
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
 *
 *******************************************************************************
 *
 * @return
 *          \retval MORSE_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa MORSE_ztpgqrt_Tile
 * @sa MORSE_ztpgqrt_Tile_Async
 * @sa MORSE_ctpgqrt
 * @sa MORSE_dtpgqrt
 * @sa MORSE_stpgqrt
 * @sa MORSE_zgeqrs
 *
 ******************************************************************************/
int MORSE_ztpgqrt( int M, int N, int K, int L,
136 137 138 139
                   MORSE_Complex64_t *V1, int LDV1, MORSE_desc_t *descT1,
                   MORSE_Complex64_t *V2, int LDV2, MORSE_desc_t *descT2,
                   MORSE_Complex64_t *Q1, int LDQ1,
                   MORSE_Complex64_t *Q2, int LDQ2 )
140 141 142 143 144 145
{
    int NB;
    int status;
    MORSE_context_t *morse;
    MORSE_sequence_t *sequence = NULL;
    MORSE_request_t request = MORSE_REQUEST_INITIALIZER;
146 147 148 149
    MORSE_desc_t descQ1l, descQ1t;
    MORSE_desc_t descQ2l, descQ2t;
    MORSE_desc_t descV1l, descV1t;
    MORSE_desc_t descV2l, descV2t;
150
    int minMK = chameleon_min( M, K );
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174

    morse = morse_context_self();
    if (morse == NULL) {
        morse_fatal_error("MORSE_ztpgqrt", "MORSE not initialized");
        return MORSE_ERR_NOT_INITIALIZED;
    }

    /* Check input arguments */
    if (M < 0) {
        morse_error("MORSE_ztpgqrt", "illegal value of M");
        return -1;
    }
    if (N < 0) {
        morse_error("MORSE_ztpgqrt", "illegal value of N");
        return -2;
    }
    if (K < 0) {
        morse_error("MORSE_ztpgqrt", "illegal value of K");
        return -3;
    }
    if ((L < 0) || ((L > minMK) && (minMK > 0))) {
        morse_error("MORSE_ztpgqrt", "illegal value of N");
        return -4;
    }
175 176 177 178 179 180
    if (K != N) {
        morse_error("MORSE_ztpgqrt", "illegal value of K and N. K must be equal to N");
        return -3;
    }
    if (LDV1 < chameleon_max(1, K)) {
        morse_error("MORSE_ztpgqrt", "illegal value of LDV1");
181 182
        return -6;
    }
183 184
    if (LDV2 < chameleon_max(1, M)) {
        morse_error("MORSE_ztpgqrt", "illegal value of LDV2");
185 186
        return -9;
    }
187 188
    if (LDQ1 < chameleon_max(1, K)) {
        morse_error("MORSE_ztpgqrt", "illegal value of LDQ1");
189 190
        return -11;
    }
191 192 193 194
    if (LDQ2 < chameleon_max(1, M)) {
        morse_error("MORSE_ztpgqrt", "illegal value of LDQ2");
        return -13;
    }
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209

    /* Quick return */
    if (minMK == 0)
        return MORSE_SUCCESS;

    /* Tune NB & IB depending on M, N & NRHS; Set NBNBSIZE */
    status = morse_tune(MORSE_FUNC_ZGELS, M, K, 0);
    if (status != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt", "morse_tune() failed");
        return status;
    }

    /* Set NT */
    NB = MORSE_NB;

210
    morse_sequence_create( morse, &sequence );
211

212
    /* Submit the matrix conversion */
213
    morse_zlap2tile( morse, &descV1l, &descV1t, MorseDescInput, MorseUpperLower,
214
                     V1, NB, NB, LDV1, K, M, K, sequence, &request );
215
    morse_zlap2tile( morse, &descV2l, &descV2t, MorseDescInput, MorseUpperLower,
216
                     V2, NB, NB, LDV2, K, M, K, sequence, &request );
217
    morse_zlap2tile( morse, &descQ1l, &descQ1t, MorseDescInout, MorseUpperLower,
218
                     Q1, NB, NB, LDQ1, N, K, N, sequence, &request );
219
    morse_zlap2tile( morse, &descQ2l, &descQ2t, MorseDescInout, MorseUpperLower,
220
                     Q2, NB, NB, LDQ2, N, M, N, sequence, &request );
221 222

    /* Call the tile interface */
223
    MORSE_ztpgqrt_Tile_Async( L, &descV1t, descT1, &descV2t, descT2, &descQ1t, &descQ2t, sequence, &request );
224

Mathieu Faverge's avatar
Mathieu Faverge committed
225
    /* Submit the matrix conversion back */
226 227 228 229
    morse_ztile2lap( morse, &descV1l, &descV1t,
                     MorseDescInput, MorseUpperLower, sequence, &request );
    morse_ztile2lap( morse, &descV2l, &descV2t,
                     MorseDescInput, MorseUpperLower, sequence, &request );
Mathieu Faverge's avatar
Mathieu Faverge committed
230
    morse_ztile2lap( morse, &descQ1l, &descQ1t,
231
                     MorseDescInout, MorseUpperLower, sequence, &request );
Mathieu Faverge's avatar
Mathieu Faverge committed
232
    morse_ztile2lap( morse, &descQ2l, &descQ2t,
233
                     MorseDescInout, MorseUpperLower, sequence, &request );
234 235
    MORSE_Desc_Flush( descT1, sequence );
    MORSE_Desc_Flush( descT2, sequence );
Mathieu Faverge's avatar
Mathieu Faverge committed
236

237
    morse_sequence_wait( morse, sequence );
Mathieu Faverge's avatar
Mathieu Faverge committed
238

Mathieu Faverge's avatar
Mathieu Faverge committed
239 240 241 242 243
    /* Cleanup the temporary data */
    morse_ztile2lap_cleanup( morse, &descV1l, &descV1t );
    morse_ztile2lap_cleanup( morse, &descV2l, &descV2t );
    morse_ztile2lap_cleanup( morse, &descQ1l, &descQ1t );
    morse_ztile2lap_cleanup( morse, &descQ2l, &descQ2t );
244 245

    status = sequence->status;
246
    morse_sequence_destroy( morse, sequence );
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
    return status;
}

/**
 *******************************************************************************
 *
 * @ingroup MORSE_Complex64_t_Tile
 *
 *  MORSE_ztpgqrt_Tile - Generates a partial Q matrix formed with a blocked QR
 *  factorization of a "triangular-pentagonal" matrix C, which is composed of an
 *  unused triangular block and a pentagonal block V, using the compact
 *  representation for Q. See MORSE_ztpqrt() to generate V.
 *
 *******************************************************************************
 *
 *
 *******************************************************************************
 *
 * @return
 *          \retval MORSE_SUCCESS successful exit
 *
 *******************************************************************************
 *
 * @sa MORSE_ztpgqrt
 * @sa MORSE_ztpgqrt_Tile_Async
 * @sa MORSE_ctpgqrt_Tile
 * @sa MORSE_dtpgqrt_Tile
 * @sa MORSE_stpgqrt_Tile
 * @sa MORSE_zgeqrs_Tile
 *
 ******************************************************************************/
278 279 280 281
int MORSE_ztpgqrt_Tile( int L,
                        MORSE_desc_t *V1, MORSE_desc_t *T1,
                        MORSE_desc_t *V2, MORSE_desc_t *T2,
                        MORSE_desc_t *Q1, MORSE_desc_t *Q2 )
282 283 284 285 286 287 288 289 290 291 292
{
    MORSE_context_t *morse;
    MORSE_sequence_t *sequence = NULL;
    MORSE_request_t request = MORSE_REQUEST_INITIALIZER;
    int status;

    morse = morse_context_self();
    if (morse == NULL) {
        morse_fatal_error("MORSE_ztpgqrt_Tile", "MORSE not initialized");
        return MORSE_ERR_NOT_INITIALIZED;
    }
293 294
    morse_sequence_create( morse, &sequence );

295
    MORSE_ztpgqrt_Tile_Async( L, V1, T1, V2, T2, Q1, Q2, sequence, &request );
296 297 298 299 300

    MORSE_Desc_Flush( V1, sequence );
    MORSE_Desc_Flush( T1, sequence );
    MORSE_Desc_Flush( V2, sequence );
    MORSE_Desc_Flush( T2, sequence );
301 302
    MORSE_Desc_Flush( Q1, sequence );
    MORSE_Desc_Flush( Q2, sequence );
Mathieu Faverge's avatar
Mathieu Faverge committed
303

304
    morse_sequence_wait( morse, sequence );
305
    status = sequence->status;
306
    morse_sequence_destroy( morse, sequence );
307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338
    return status;
}

/**
 *******************************************************************************
 *
 * @ingroup MORSE_Complex64_t_Tile_Async
 *
 *  MORSE_ztpgqrt_Tile_Async - Generates a partial Q matrix formed with a blocked QR
 *  factorization of a "triangular-pentagonal" matrix C, which is composed of an
 *  unused triangular block and a pentagonal block V, using the compact
 *  representation for Q. See MORSE_ztpqrt() to generate V.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa MORSE_ztpgqrt
 * @sa MORSE_ztpgqrt_Tile
 * @sa MORSE_ctpgqrt_Tile_Async
 * @sa MORSE_dtpgqrt_Tile_Async
 * @sa MORSE_stpgqrt_Tile_Async
 * @sa MORSE_zgeqrs_Tile_Async
 *
 ******************************************************************************/
339 340 341 342
int MORSE_ztpgqrt_Tile_Async( int L,
                              MORSE_desc_t *V1, MORSE_desc_t *T1,
                              MORSE_desc_t *V2, MORSE_desc_t *T2,
                              MORSE_desc_t *Q1, MORSE_desc_t *Q2,
343 344 345
                              MORSE_sequence_t *sequence, MORSE_request_t *request )
{
    MORSE_context_t *morse;
346
    MORSE_desc_t D, *Dptr = NULL;
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361

    morse = morse_context_self();
    if (morse == NULL) {
        morse_error("MORSE_ztpgqrt_Tile", "MORSE not initialized");
        return MORSE_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        morse_fatal_error("MORSE_ztpgqrt_Tile", "NULL sequence");
        return MORSE_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        morse_fatal_error("MORSE_ztpgqrt_Tile", "NULL request");
        return MORSE_ERR_UNALLOCATED;
    }
    /* Check sequence status */
Mathieu Faverge's avatar
Mathieu Faverge committed
362
    if (sequence->status == MORSE_SUCCESS) {
363
        request->status = MORSE_SUCCESS;
Mathieu Faverge's avatar
Mathieu Faverge committed
364 365
    }
    else {
366
        return morse_request_fail(sequence, request, MORSE_ERR_SEQUENCE_FLUSHED);
Mathieu Faverge's avatar
Mathieu Faverge committed
367
    }
368 369

    /* Check descriptors for correctness */
370 371 372 373 374 375 376 377 378 379
    if (morse_desc_check(V1) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid V1 descriptor");
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
    if (morse_desc_check(T1) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid T1 descriptor");
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
    if (morse_desc_check(V2) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid V2 descriptor");
380 381
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
382 383
    if (morse_desc_check(T2) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid T2 descriptor");
384 385
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
386 387
    if (morse_desc_check(Q1) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid Q1 descriptor");
388 389
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
390 391
    if (morse_desc_check(Q2) != MORSE_SUCCESS) {
        morse_error("MORSE_ztpgqrt_Tile", "invalid Q2 descriptor");
392 393 394
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
    /* Check input arguments */
395
    if (Q1->nb != Q1->mb) {
396 397 398
        morse_error("MORSE_ztpgqrt_Tile", "only square tiles supported");
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
399
    if ( (L != 0) && (((Q2->m - L) % Q2->mb) != 0) ) {
400 401 402
        morse_error("MORSE_ztpgqrt_Tile", "Triangular part must be aligned with tiles");
        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
    }
403 404 405 406
#if defined(CHAMELEON_COPY_DIAG)
    {
        int minMT;
        if (V1->m > V1->n) {
407
            minMT = V1->nt;
408 409 410 411 412 413 414
        } else {
            minMT = V1->mt;
        }
        morse_zdesc_alloc_diag(D, V1->mb, V1->nb, minMT*V1->mb, V1->nb, 0, 0, minMT*V1->mb, V1->nb, V1->p, V1->q);
        Dptr = &D;
    }
#endif
415 416

    /* if (morse->householder == MORSE_FLAT_HOUSEHOLDER) { */
417 418
    morse_pzlaset( MorseUpperLower, 0., 1., Q1, sequence, request );
    morse_pzlaset( MorseUpperLower, 0., 0., Q2, sequence, request );
419
    morse_pztpgqrt( L, V1, T1, V2, T2, Q1, Q2, Dptr, sequence, request );
420 421
    /* } */
    /* else { */
Mathieu Faverge's avatar
Mathieu Faverge committed
422
    /*    morse_pztpgqrtrh( Q1, T, MORSE_RHBLK, sequence, request ); */
423
    /* } */
424
    if (Dptr != NULL) {
425 426 427 428 429 430 431 432
        MORSE_Desc_Flush( V1, sequence );
        MORSE_Desc_Flush( T1, sequence );
        MORSE_Desc_Flush( V2, sequence );
        MORSE_Desc_Flush( T2, sequence );
        MORSE_Desc_Flush( Q1, sequence );
        MORSE_Desc_Flush( Q2, sequence );
        MORSE_Desc_Flush( Dptr, sequence );
        morse_sequence_wait( morse, sequence );
Mathieu Faverge's avatar
Mathieu Faverge committed
433
        morse_desc_mat_free( Dptr );
434
    }
Mathieu Faverge's avatar
Mathieu Faverge committed
435
    (void)D;
436 437
    return MORSE_SUCCESS;
}