My professional and academic background is in numerical analysis and scientific computing (including BLAS/LAPACK level implementations), but I admit I haven't done deep numerical implementations in years.
At the other end of the spectrum, getting a matrix-matrix multiply to run fast isn’t easy either. It’s what necessitated the kind of approach the authors of BLIS adopted. On paper it’s easy, but actually getting it to run fast on a computer isn’t.
And then, there are many parts of the algorithms deep in things like the various matrix decompositions that are difficult to parallelise because of non-trivial data dependencies. It’s easy to write some code to pivot a matrix; it’s very hard to do it in an efficient and scalable way. And because all the higher-level routines depend on them, so they have a large effect on overall performance.
With that out of the way, basic linear algebra operations that require sophisticated algorithms and are not "embarrassingly parallel": matrix multiplication, matrix inversion, matrix decomposition (SVM, QR, etc). Some of these algos fall into BLAS and others in LAPACK.
As soon as you start mixing numerical stability (data-dependant operation reordering) concerns into linear algebra routines -- arguably one of the main points of LAPACK -- parallelization necessarily gets trickier because you either need a good enough heuristic to ignore those reorderings or you also need data-dependant scheduling without too many pathologies.
But from this exercise alone, knowing the tricks you used already, you can see how un-embarrassingly parallel this task is (frankly if it is truly embarrassingly parallel, then `#prama omp for` should bring you close to best possible performance already.)
I don’t think your prof got it wrong, it is just that you misunderstood them.