Merge back qdwh to integrate lascal

compute/CMakeLists.txt

@@ -66,8 +66,10 @@ set(ZSRC
     pztrsm.c
     pztrsmpl.c
     pztradd.c
+    pzlascal.c
     ###
     zgeadd.c
+    zlascal.c
     zgemm.c
     zhemm.c
     zher2k.c

compute/pzlascal.c

+/**
+ *
+ * @file pzlascal.c
+ *
+ *  MORSE auxiliary routines
+ *  MORSE is a software package provided by Univ. of Tennessee,
+ *  Univ. of California Berkeley and Univ. of Colorado Denver
+ *
+ * @version 2.8.0
+ * @author Dalal Sukkari
+ * @date 2010-11-15
+ * @precisions normal z -> s d c
+ *
+ **/
+#include "control/common.h"
+
+#define A(m, n) A,  m,  n
+/**
+ *  Parallel scale of a matrix A
+ **/
+void morse_pzlascal(MORSE_enum uplo, MORSE_Complex64_t alpha, MORSE_desc_t *A,
+                    MORSE_sequence_t *sequence, MORSE_request_t *request)
+{
+    MORSE_context_t *morse;
+    MORSE_option_t options;
+
+    int tempmm, tempnn, tempmn, tempnm;
+    int m, n;
+    int ldam, ldan;
+    int minmnt = min(A->mt, A->nt);
+
+    morse = morse_context_self();
+    if (sequence->status != MORSE_SUCCESS)
+        return;
+
+    RUNTIME_options_init(&options, morse, sequence, request);
+
+    switch(uplo) {
+    case MorseLower:
+        for (n = 0; n < minmnt; n++) {
+            tempnm = n == A->mt-1 ? A->m-n*A->mb : A->mb;
+            tempnn = n == A->nt-1 ? A->n-n*A->nb : A->nb;
+            ldan = BLKLDD(A, n);
+
+            MORSE_TASK_zlascal(
+                &options,
+                MorseLower, tempnm, tempnn, A->mb,
+                alpha, A(n, n), ldan);
+
+            for (m = n+1; m < A->mt; m++) {
+                tempmm = m == A->mt-1 ? A->m-A->mb*m : A->nb;
+                ldam = BLKLDD(A, m);
+
+                MORSE_TASK_zlascal(
+                    &options,
+                    MorseUpperLower, tempmm, tempnn, A->mb,
+                    alpha, A(m, n), ldam);
+            }
+        }
+        break;
+
+    case MorseUpper:
+        for (m = 0; m < minmnt; m++) {
+            tempmm = m == A->mt-1 ? A->m-A->mb*m : A->nb;
+            tempmn = m == A->nt-1 ? A->n-m*A->nb : A->nb;
+            ldam = BLKLDD(A, m);
+
+            MORSE_TASK_zlascal(
+                &options,
+                MorseUpper, tempmm, tempmn, A->mb,
+                alpha, A(m, m), ldam);
+
+            for (n = m+1; n < A->nt; n++) {
+                tempnn = n == A->nt-1 ? A->n-n*A->nb : A->nb;
+
+                MORSE_TASK_zlascal(
+                    &options,
+                    MorseUpperLower, tempmm, tempnn, A->mb,
+                    alpha, A(m, n), ldam);
+            }
+        }
+        break;
+
+    case MorseUpperLower:
+    default:
+        for (m = 0; m < A->mt; m++) {
+            tempmm = m == A->mt-1 ? A->m-A->mb*m : A->nb;
+            ldam = BLKLDD(A, m);
+
+            for (n = 0; n < A->nt; n++) {
+                tempnn = n == A->nt-1 ? A->n-n*A->nb : A->nb;
+
+                MORSE_TASK_zlascal(
+                    &options,
+                    MorseUpperLower, tempmm, tempnn, A->mb,
+                    alpha, A(m, n), ldam);
+            }
+        }
+    }
+    RUNTIME_options_ws_free(&options);
+    RUNTIME_options_finalize(&options, morse);
+    MORSE_TASK_dataflush_all();
+}

compute/zlascal.c

+/**
+ *
+ * @file zlascal.c
+ *
+ *  MORSE computational routines
+ *  MORSE is a software package provided by Univ. of Tennessee,
+ *  Univ. of California Berkeley and Univ. of Colorado Denver
+ *
+ * @version 2.8.0
+ * @author Dalal Sukkari
+ * @date 2010-11-15
+ * @precisions normal z -> s d c
+ *
+ **/
+#include "control/common.h"
+
+/***************************************************************************//**
+ *
+ * @ingroup MORSE_Complex64_t
+ *
+ *  MORSE_zlascal - Scales a matrix by the scalar alpha as in
+ *  ScaLAPACK pzlascal().
+ *
+ *    \f[ A = \alpha A \f],
+ *
+ *  alpha is a scalar, and A a general, upper or lower trapezoidal matrix.
+ *
+ *******************************************************************************
+ *
+ * @param[in] uplo
+ *          Specifies the shape of A:
+ *          = MorseUpperLower: A is a general matrix.
+ *          = MorseUpper: A is an upper trapezoidal matrix.
+ *          = MorseLower: A is a lower trapezoidal matrix.
+ *
+ * @param[in] M
+ *          M specifies the number of rows of the matrix A. M >= 0.
+ *
+ * @param[in] N
+ *          N specifies the number of columns of the matrix A. N >= 0.
+ *
+ * @param[in] alpha
+ *          alpha specifies the scalar alpha
+ *
+ * @param[in,out] A
+ *          A is a LDA-by-N matrix.
+ *
+ * @param[in] LDA
+ *          The leading dimension of the array A. LDA >= max(1,M).
+ *
+ *******************************************************************************
+ *
+ * @return
+ *          \retval MORSE_SUCCESS successful exit
+ *
+ *******************************************************************************
+ *
+ * @sa MORSE_zlascal_Tile
+ * @sa MORSE_clascal
+ * @sa MORSE_dlascal
+ * @sa MORSE_slascal
+ *
+ ******************************************************************************/
+int MORSE_zlascal(MORSE_enum uplo, int M, int N,
+                   MORSE_Complex64_t alpha, MORSE_Complex64_t *A, int LDA)
+{
+    int NB;
+    int status;
+    MORSE_desc_t descA;
+    MORSE_context_t *morse;
+    MORSE_sequence_t *sequence = NULL;
+    MORSE_request_t request = MORSE_REQUEST_INITIALIZER;
+
+    morse = morse_context_self();
+    if (morse == NULL) {
+        morse_fatal_error("MORSE_zlascal", "MORSE not initialized");
+        return MORSE_ERR_NOT_INITIALIZED;
+    }
+
+    /* Check input arguments */
+    if (uplo != MorseUpper && uplo != MorseLower && uplo != MorseUpperLower) {
+        morse_error("MORSE_zlascal", "illegal value of uplo");
+        return -1;
+    }
+    if (M < 0) {
+        morse_error("MORSE_zlascal", "illegal value of M");
+        return -2;
+    }
+    if (N < 0) {
+        morse_error("MORSE_zlascal", "illegal value of N");
+        return -3;
+    }
+    if (LDA < max(1, M)) {
+        morse_error("MORSE_zlascal", "illegal value of LDA");
+        return -6;
+    }
+
+    /* Quick return */
+    if (M == 0 || N == 0 ||
+        (alpha == (MORSE_Complex64_t)1.0))
+        return MORSE_SUCCESS;
+
+    /* Tune NB depending on M, N & NRHS; Set NBNBSIZE */
+    status = morse_tune(MORSE_FUNC_ZGEMM, M, N, 0);
+    if (status != MORSE_SUCCESS) {
+        morse_error("MORSE_zlascal", "morse_tune() failed");
+        return status;
+    }
+
+    /* Set MT & NT & KT */
+    NB = MORSE_NB;
+
+    morse_sequence_create(morse, &sequence);
+
+/*    if ( MORSE_TRANSLATION == MORSE_OUTOFPLACE ) {*/
+        morse_zooplap2tile( descA, A, NB, NB, LDA, N, 0, 0, M, N, sequence, &request,
+                             morse_desc_mat_free(&(descA)) );
+/*    } else {*/
+/*        morse_ziplap2tile( descA, A, NB, NB, LDA, N , 0, 0, M, N,*/
+/*                            sequence, &request);*/
+/*    }*/
+
+    /* Call the tile interface */
+    MORSE_zlascal_Tile_Async(
+        uplo, alpha, &descA, sequence, &request);
+
+/*    if ( MORSE_TRANSLATION == MORSE_OUTOFPLACE ) {*/
+        morse_zooptile2lap( descA, A, NB, NB, LDA, N,  sequence, &request);
+        RUNTIME_sequence_wait(morse, sequence);
+        morse_desc_mat_free(&descA);
+/*    } else {*/
+/*        morse_ziptile2lap( descA, A, NB, NB, LDA, N, sequence, &request);*/
+/*        morse_dynamic_sync();*/
+/*    }*/
+
+    status = sequence->status;
+    morse_sequence_destroy(morse, sequence);
+    return status;
+}
+
+/***************************************************************************//**
+ *
+ * @ingroup MORSE_Complex64_t_Tile
+ *
+ *  MORSE_zlascal_Tile - Scales a matrix by the scalar alpha as in
+ *  ScaLAPACK pzlascal().
+ *
+ *    \f[ A = \alpha A \f],
+ *
+ *  alpha is a scalar, and A a general, upper or lower trapezoidal matrix.
+ *
+ *******************************************************************************
+ *
+ * @param[in] uplo
+ *          Specifies the shape of A:
+ *          = MorseUpperLower: A is a general matrix.
+ *          = MorseUpper: A is an upper trapezoidal matrix.
+ *          = MorseLower: A is a lower trapezoidal matrix.
+ *
+ * @param[in] alpha
+ *          alpha specifies the scalar alpha
+ *
+ * @param[in] A
+ *          A is a LDA-by-N matrix.
+ *
+ *******************************************************************************
+ *
+ * @return
+ *          \retval MORSE_SUCCESS successful exit
+ *
+ *******************************************************************************
+ *
+ * @sa MORSE_zlascal
+ * @sa MORSE_zlascal_Tile_Async
+ * @sa MORSE_clascal_Tile
+ * @sa MORSE_dlascal_Tile
+ * @sa MORSE_slascal_Tile
+ *
+ ******************************************************************************/
+int MORSE_zlascal_Tile(MORSE_enum uplo,
+                        MORSE_Complex64_t alpha, MORSE_desc_t *A)
+{
+    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_zlascal_Tile", "MORSE not initialized");
+        return MORSE_ERR_NOT_INITIALIZED;
+    }
+    morse_sequence_create(morse, &sequence);
+    MORSE_zlascal_Tile_Async(uplo, alpha, A, sequence, &request);
+    RUNTIME_sequence_wait(morse, sequence);
+    status = sequence->status;
+    morse_sequence_destroy(morse, sequence);
+    return status;
+}
+
+/***************************************************************************//**
+ *
+ * @ingroup MORSE_Complex64_t_Tile_Async
+ *
+ *  MORSE_zlascal_Tile_Async - Scales a matrix by the scalar alpha as in
+ *  ScaLAPACK pzlascal().
+ *  Non-blocking equivalent of MORSE_zlascal_Tile().
+ *  May return before the computation is finished.
+ *  Allows for pipelining of operations at runtime.
+ *
+ *******************************************************************************
+ *
+ * @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_zlascal
+ * @sa MORSE_zlascal_Tile
+ * @sa MORSE_clascal_Tile_Async
+ * @sa MORSE_dlascal_Tile_Async
+ * @sa MORSE_slascal_Tile_Async
+ *
+ ******************************************************************************/
+int MORSE_zlascal_Tile_Async(MORSE_enum uplo,
+                              MORSE_Complex64_t alpha, MORSE_desc_t *A,
+                              MORSE_sequence_t *sequence, MORSE_request_t *request)
+{
+    MORSE_context_t *morse;
+    MORSE_desc_t descA;
+
+    morse = morse_context_self();
+    if (morse == NULL) {
+        morse_fatal_error("MORSE_zlascal_Tile_Async", "MORSE not initialized");
+        return MORSE_ERR_NOT_INITIALIZED;
+    }
+    if (sequence == NULL) {
+        morse_fatal_error("MORSE_zlascal_Tile_Async", "NULL sequence");
+        return MORSE_ERR_UNALLOCATED;
+    }
+    if (request == NULL) {
+        morse_fatal_error("MORSE_zlascal_Tile_Async", "NULL request");
+        return MORSE_ERR_UNALLOCATED;
+    }
+    /* Check sequence status */
+    if (sequence->status == MORSE_SUCCESS)
+        request->status = MORSE_SUCCESS;
+    else
+        return morse_request_fail(sequence, request, MORSE_ERR_SEQUENCE_FLUSHED);
+
+    /* Check descriptors for correctness */
+    if (morse_desc_check(A) != MORSE_SUCCESS) {
+        morse_error("MORSE_zlascal_Tile_Async", "invalid first descriptor");
+        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
+    } else {
+        descA = *A;
+    }
+    /* Check input arguments */
+    if (uplo != MorseUpper && uplo != MorseLower && uplo != MorseUpperLower) {
+        morse_error("MORSE_zlascal", "illegal value of uplo");
+        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
+    }
+
+    if ( (descA.i%descA.mb != 0) || (descA.j%descA.nb != 0) ) {
+        morse_error("MORSE_zlascal", "start indexes have to be multiple of tile size");
+        return morse_request_fail(sequence, request, MORSE_ERR_ILLEGAL_VALUE);
+    }
+
+    /* Quick return */
+    if ( (descA.m == 0) || (descA.n == 0) ||
+         (alpha == (MORSE_Complex64_t)1.0) )
+        return MORSE_SUCCESS;
+
+    morse_pzlascal( uplo, alpha, A, sequence, request);
+
+    return MORSE_SUCCESS;
+}

control/compute_z.h

@@ -118,6 +118,7 @@ void morse_pzlanhe(MORSE_enum norm, MORSE_enum uplo, MORSE_desc_t *A, double *re
 #endif
 void morse_pzlansy(MORSE_enum norm, MORSE_enum uplo, MORSE_desc_t *A, double *result, MORSE_sequence_t *sequence, MORSE_request_t *request);
 void morse_pzlantr(MORSE_enum norm, MORSE_enum uplo, MORSE_enum diag, MORSE_desc_t *A, double *result, MORSE_sequence_t *sequence, MORSE_request_t *request);
+void morse_pzlascal(MORSE_enum uplo, MORSE_Complex64_t alpha, MORSE_desc_t *A, MORSE_sequence_t *sequence, MORSE_request_t *request);
 void morse_pzlaset( MORSE_enum uplo, MORSE_Complex64_t alpha, MORSE_Complex64_t beta, MORSE_desc_t *A, MORSE_sequence_t *sequence, MORSE_request_t *request);
 void morse_pzlaset2(MORSE_enum uplo, MORSE_Complex64_t alpha,                          MORSE_desc_t *A, MORSE_sequence_t *sequence, MORSE_request_t *request);
 void morse_pzlaswp(MORSE_desc_t *B, int *IPIV, int inc, MORSE_sequence_t *sequence, MORSE_request_t *request);

coreblas/compute/CMakeLists.txt

@@ -33,6 +33,7 @@ set(ZSRC
     core_dzasum.c
     core_zaxpy.c
     core_zgeadd.c
+    core_zlascal.c
     core_zgelqt.c
     core_zgemm.c
     core_zgeqrt.c

coreblas/compute/core_zlascal.c

+/**
+ * @file core_zlascal.c
+ *
+ *  MORSE computational routines
+ *  MORSE is a software package provided by Univ. of Tennessee,
+ *  Univ. of California Berkeley and Univ. of Colorado Denver
+ *
+ * @version 2.8.0
+ * @author Dalal Sukkari
+ * @date 2015-11-05
+ * @precisions normal z -> c d s
+ *
+ **/
+#include "coreblas/include/coreblas.h"
+#include "coreblas/include/cblas.h"
+#include <math.h>
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup CORE_MORSE_Complex64_t
+ *
+ *  CORE_zlascal scales a two-dimensional matrix A. As opposite to
+ *  CORE_zlascl(), no checks is performed to prevent under/overflow. This should
+ *  have been done at higher level.
+ *
+ *******************************************************************************
+ *
+ * @param[in] uplo
+ *          Specifies the shape of A:
+ *          = MorseUpperLower: A is a general matrix.
+ *          = MorseUpper: A is an upper trapezoidal matrix.
+ *          = MorseLower: A is a lower trapezoidal matrix.
+ *
+ * @param[in] m is the number of rows of the matrix A. m >= 0
+ *
+ * @param[in] n is the number of columns of the matrix A. n >= 0
+ *
+ * @param[in] alpha
+ *            The scalar factor.
+ *
+ * @param[in,out] A is the matrix to be multiplied by alpha
+ *
+ * @param[in] lda is the leading dimension of the array A. lda >= max(1,m).
+ *
+ *******************************************************************************
+ *
+ * @return
+ *          \retval MORSE_SUCCESS successful exit
+ *          \retval <0 if -i, the i-th argument had an illegal value
+ *
+ ******************************************************************************/
+int
+CORE_zlascal( MORSE_enum uplo, int m, int n,
+              MORSE_Complex64_t alpha, MORSE_Complex64_t *A, int lda )
+{
+    int i;
+
+    if ( (uplo != MorseUpperLower) &&
+         (uplo != MorseUpper)      &&
+         (uplo != MorseLower))
+    {
+        coreblas_error(1, "illegal value of uplo");
+        return -1;
+    }
+
+    if (m < 0) {
+        coreblas_error(2, "Illegal value of m");
+        return -2;
+    }
+    if (n < 0) {
+        coreblas_error(3, "Illegal value of n");
+        return -3;
+    }
+    if ( (lda < max(1,m)) && (m > 0) ) {
+        coreblas_error(6, "Illegal value of lda");
+        return -6;
+    }
+
+    switch ( uplo ) {
+    case MorseUpper:
+        for(i=0; i<n; i++) {
+            cblas_zscal( min( i+1, m ), CBLAS_SADDR(alpha), A+i*lda, 1 );
+        }
+        break;
+
+    case MorseLower:
+        for(i=0; i<n; i++) {
+            cblas_zscal( max( m, m-i ), CBLAS_SADDR(alpha), A+i*lda, 1 );
+        }
+        break;
+    default:
+        if (m == lda) {
+            cblas_zscal( m*n, CBLAS_SADDR(alpha), A, 1 );
+        }
+        else {
+            for(i=0; i<n; i++) {
+                cblas_zscal( m, CBLAS_SADDR(alpha), A+i*lda, 1 );
+            }
+        }
+    }
+
+    return MORSE_SUCCESS;
+}

coreblas/compute/core_ztsmlq.c

@@ -259,10 +259,10 @@ int CORE_ztsmlq(MORSE_enum side, MORSE_enum trans,
         CORE_zparfb(
             side, trans, MorseForward, MorseRowwise,
             mi, ni, M2, N2, kb, 0,
-            &A1[LDA1*jc+ic], LDA1,
+            A1 + LDA1 * jc + ic, LDA1,
             A2, LDA2,
-            &V[i], LDV,
-            &T[LDT*i], LDT,
+            V + i,       LDV,
+            T + i * LDT, LDT,
             WORK, LDWORK);
     }
     return MORSE_SUCCESS;

coreblas/include/coreblas_z.h

@@ -63,6 +63,8 @@ int CORE_zgeadd(MORSE_enum trans, int M, int N,
                 const MORSE_Complex64_t *A, int LDA,
                       MORSE_Complex64_t beta,
                       MORSE_Complex64_t *B, int LDB);
+int CORE_zlascal( MORSE_enum uplo, int m, int n,
+              MORSE_Complex64_t alpha, MORSE_Complex64_t *A, int lda );
 int  CORE_zgelqt(int M, int N, int IB,
                  MORSE_Complex64_t *A, int LDA,
                  MORSE_Complex64_t *T, int LDT,

cudablas/compute/cuda_zgelqt.c

@@ -26,17 +26,17 @@
 
 #if defined(CHAMELEON_USE_MAGMA)
 int CUDA_zgelqt(
-        magma_int_t m, magma_int_t n, magma_int_t nb,
-        magmaDoubleComplex *da, magma_int_t ldda,
-        magmaDoubleComplex *v,  magma_int_t ldv,
-        magmaDoubleComplex *dt, magma_int_t lddt,
-        magmaDoubleComplex *t,  magma_int_t ldt,
-        magmaDoubleComplex *dd,
-        magmaDoubleComplex *d,  magma_int_t ldd,
-        magmaDoubleComplex *tau,
-        magmaDoubleComplex *hwork,
-        magmaDoubleComplex *dwork,
-        CUstream stream)
+    magma_int_t m, magma_int_t n, magma_int_t nb,
+    magmaDoubleComplex *da, magma_int_t ldda,
+    magmaDoubleComplex *v,  magma_int_t ldv,
+    magmaDoubleComplex *dt, magma_int_t lddt,
+    magmaDoubleComplex *t,  magma_int_t ldt,
+    magmaDoubleComplex *dd,
+    magmaDoubleComplex *d,  magma_int_t ldd,
+    magmaDoubleComplex *tau,
+    magmaDoubleComplex *hwork,
+    magmaDoubleComplex *dwork,
+    CUstream stream)
 {
 #define da_ref(a_1,a_2) ( da+(a_2)*(ldda) + (a_1))
 #define v_ref(a_1,a_2)  ( v+(a_2)*(ldv) + (a_1))
@@ -47,17 +47,17 @@ int CUDA_zgelqt(
     double _Complex one=1.;
 
     if (m < 0) {
-    return -1;
+        return -1;
     } else if (n < 0) {
-    return -2;
+        return -2;
     } else if (ldda < max(1,m)) {
-    return -4;
+        return -4;
     }
 
     k = min(m,n);
     if (k == 0) {
-    hwork[0] = *((magmaDoubleComplex*) &one);
-    return MAGMA_SUCCESS;
+        hwork[0] = *((magmaDoubleComplex*) &one);
+        return MAGMA_SUCCESS;
     }
 
     /* lower parts of little T must be zero: memset to 0 for simplicity */
@@ -92,7 +92,7 @@ int CUDA_zgelqt(
 
             magma_queue_sync( stream );
             /* Form the triangular factor of the block reflector on the host
-            H = H'(i+ib-1) . . . H(i+1) H(i) */
+             H = H'(i+ib-1) . . . H(i+1) H(i) */
             CORE_zgelqt(ib, cols, ib,
                         (double _Complex*) v_ref(0,i), ib,
                         (double _Complex*) t_ref(0,0), ib,
@@ -100,7 +100,7 @@ int CUDA_zgelqt(
                         (double _Complex*) hwork);
 
             /* put 0s in the lower triangular part of a panel (and 1s on the
-              diagonal); copy the lower triangular in d */
+             diagonal); copy the lower triangular in d */
             CORE_zgesplit(MorseRight, MorseUnit, ib, min(ib,cols),
                           (double _Complex*) v_ref(0,i), ib,
                           (double _Complex*) d, ib);