pzgebrd_ge2gb.c 4.45 KB
Newer Older
1
/**
2 3
 *
 * @file pzgebrd_ge2gb.c
4
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
5 6
 * @copyright 2009-2014 The University of Tennessee and The University of
 *                      Tennessee Research Foundation. All rights reserved.
7 8
 * @copyright 2012-2014 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
9
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
10
 ***
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
 *
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Hatem Ltaief
 * @author Azzam Haidar
 * @date 2010-11-15
 * @precisions normal z -> s d c
 *
 **/
#include "control/common.h"

BOUCHERIE Raphael's avatar
BOUCHERIE Raphael committed
26
void morse_pzgebrd_ge2gb(MORSE_desc_t A, MORSE_desc_t T, MORSE_desc_t D,
27 28 29 30 31 32 33 34 35 36 37 38
                         MORSE_sequence_t *sequence, MORSE_request_t *request)
{
    int k;
    int tempkm, tempkn;
    if (A.m >= A.n){
       for (k = 0; k < A.nt; k++) {
           tempkm = k == A.mt-1 ? A.m-k*A.mb : A.mb;
           tempkn = k == A.nt-1 ? A.n-k*A.nb : A.nb;

           morse_pzgeqrf(
               morse_desc_submatrix(&A, k*A.mb, k*A.nb, A.m-k*A.mb, tempkn),
               morse_desc_submatrix(&T, k*T.mb, k*T.nb, T.m-k*T.mb, tempkn),
39
               morse_desc_submatrix(&D, k*T.mb, k*T.nb, T.m-k*T.mb, tempkn),
40 41 42 43 44 45 46 47
               sequence, request);

           morse_pzunmqr(
               MorseLeft,
               MorseConjTrans,
               morse_desc_submatrix(&A, k*A.mb,     k*A.nb, A.m-k*A.mb, tempkn),
               morse_desc_submatrix(&A, k*A.mb, (k+1)*A.nb, A.m-k*A.mb, A.n-(k+1)*A.nb),
               morse_desc_submatrix(&T, k*T.mb,     k*T.nb, T.m-k*T.mb, tempkn),
48
               morse_desc_submatrix(&D, k*T.mb,     k*T.nb, T.m-k*T.mb, tempkn),
49 50 51 52 53 54 55 56
               sequence, request);

           if (k+1 < A.nt){
              tempkn = k+1 == A.nt-1 ? A.n-(k+1)*A.nb : A.nb;

              morse_pzgelqf(
                  morse_desc_submatrix(&A, k*A.mb, (k+1)*A.nb, tempkm, A.n-(k+1)*A.nb),
                  morse_desc_submatrix(&T, k*T.mb, (k+1)*T.nb, T.mb,   T.n-(k+1)*T.nb),
57
                  morse_desc_submatrix(&D, k*T.mb, (k+1)*T.nb, T.mb,   T.n-(k+1)*T.nb),
58 59 60 61 62 63 64
                  sequence, request);

              morse_pzunmlq(
                  MorseRight, MorseConjTrans,
                  morse_desc_submatrix(&A,     k*A.mb, (k+1)*A.nb, tempkm,         A.n-(k+1)*A.nb),
                  morse_desc_submatrix(&A, (k+1)*A.mb, (k+1)*A.nb, A.m-(k+1)*A.mb, A.n-(k+1)*A.nb),
                  morse_desc_submatrix(&T,     k*T.mb, (k+1)*T.nb, T.mb,           T.n-(k+1)*T.nb),
65
                  morse_desc_submatrix(&D,     k*T.mb, (k+1)*T.nb, T.mb,           T.n-(k+1)*T.nb),
66 67 68 69 70 71 72 73 74 75 76 77
                  sequence, request);
           }
       }
    }
    else{
       for (k = 0; k < A.mt; k++) {
           tempkm = k == A.mt-1 ? A.m-k*A.mb : A.mb;
           tempkn = k == A.nt-1 ? A.n-k*A.nb : A.nb;

           morse_pzgelqf(
               morse_desc_submatrix(&A, k*A.mb, k*A.nb, tempkm, A.n-k*A.nb),
               morse_desc_submatrix(&T, k*T.mb, k*T.nb, T.mb,   T.n-k*T.nb),
78
               morse_desc_submatrix(&D, k*T.mb, k*T.nb, T.mb,   T.n-k*T.nb),
79 80 81 82 83 84 85
               sequence, request);

           morse_pzunmlq(
               MorseRight, MorseConjTrans,
               morse_desc_submatrix(&A,     k*A.mb, k*A.nb, tempkm,         A.n-k*A.nb),
               morse_desc_submatrix(&A, (k+1)*A.mb, k*A.nb, A.m-(k+1)*A.mb, A.n-k*A.nb),
               morse_desc_submatrix(&T,     k*T.mb, k*T.nb, T.mb,           T.n-k*T.nb),
86
               morse_desc_submatrix(&D,     k*T.mb, k*T.nb, T.mb,           T.n-k*T.nb),
87 88 89 90 91 92 93 94
               sequence, request);

           if (k+1 < A.mt){
              tempkm = k+1 == A.mt-1 ? A.m-(k+1)*A.mb : A.mb;

              morse_pzgeqrf(
                   morse_desc_submatrix(&A, (k+1)*A.mb, k*A.nb, A.m-(k+1)*A.mb, tempkn),
                   morse_desc_submatrix(&T, (k+1)*T.mb, k*T.nb, T.m-(k+1)*T.mb, tempkn),
95
                   morse_desc_submatrix(&D, (k+1)*T.mb, k*T.nb, T.m-(k+1)*T.mb, tempkn),
96 97 98 99 100 101 102
                   sequence, request);

              morse_pzunmqr(
                  MorseLeft, MorseConjTrans,
                  morse_desc_submatrix(&A, (k+1)*A.mb,     k*A.nb, A.m-(k+1)*A.mb, tempkn),
                  morse_desc_submatrix(&A, (k+1)*A.mb, (k+1)*A.nb, A.m-(k+1)*A.mb, A.n-(k+1)*A.nb),
                  morse_desc_submatrix(&T, (k+1)*T.mb,     k*T.nb, T.m-(k+1)*T.mb, tempkn),
103
                  morse_desc_submatrix(&D, (k+1)*T.mb,     k*T.nb, T.m-(k+1)*T.mb, tempkn),
104 105 106 107 108
                  sequence, request);
           }
       }
    }
}